Theory Prime_Number_Theorem.Mertens_Theorems

(*
  File:     Mertens_Theorems.thy
  Author:   Manuel Eberl (TU München)

*)
section ‹Mertens' Theorems›
theory Mertens_Theorems
imports
  Prime_Counting_Functions
  "Stirling_Formula.Stirling_Formula"
begin

(*<*)
unbundle prime_counting_notation
(*>*)

text ‹
  In this section, we will prove Mertens' First and Second Theorem. These are weaker results
  than the Prime Number Theorem, and we will derive them without using it.

  However, like Mertens himself, we will not only prove them ‹asymptotically›, but ‹absolutely›.
  This means that we will show that the remainder terms are not only ``Big-O'' of some bound, but
  we will give concrete (and reasonably tight) upper and lower bounds for them that hold on the 
  entire domain. This makes the proofs a bit more tedious.
›

subsection ‹Absolute Bounds for Mertens' First Theorem›

text ‹
  We have already shown the asymptotic form of Mertens' first theorem, i.\,e.\ 
  $\mathfrak{M}(n) = \ln n + O(1)$. We now want to obtain some absolute bounds on the $O(1)$
  remainder term using a more careful derivation than before.

  The precise bounds we will show are
  $\mathfrak {M}(n) - \ln n \in (-1-\frac{9}{\pi^2}; \ln 4] \approx (-1.9119; 1.3863]$ for
  n ∈ ℕ›.

  First, we need a simple lemma on the finiteness of exponents to consider in a sum
  of all prime powers up to a certain point:
›
lemma exponents_le_finite:
  assumes "p > (1 :: nat)" "k > 0"
  shows   "finite {i. real (p ^ (k * i + l))  x}"
proof (rule finite_subset)
  show "{i. real (p ^ (k * i + l))  x}  {..nat x}"
  proof safe
    fix i assume i: "real (p ^ (k * i + l))  x"
    have "i < 2 ^ i" by (rule less_exp)
    also from assms have "i  k * i + l" by (cases k) auto
    hence "2 ^ i  (2 ^ (k * i + l) :: nat)"
      using assms by (intro power_increasing) auto
    also have "  p ^ (k * i + l)" using assms by (intro power_mono) auto
    also have "real   x" using i by simp
    finally show "i  nat x" by linarith
  qed
qed auto

text ‹
  Next, we need the following bound on $\zeta'(2)$:
›
lemma deriv_zeta_2_bound: "Re (deriv zeta 2) > -1"
proof -
  have "((λx::real. ln (x + 3) * (x + 3) powr -2) has_integral (ln 3 + 1) / 3) (interior {0..})"
    using ln_powr_has_integral_at_top[of 1 0 3 "-2"]
    by (simp add: interior_real_atLeast powr_minus)
  hence "((λx::real. ln (x + 3) * (x + 3) powr -2) has_integral (ln 3 + 1) / 3) {0..}"
    by (subst (asm) has_integral_interior) auto
  also have "?this  ((λx::real. ln (x + 3) / (x + 3) ^ 2) has_integral (ln 3 + 1) / 3) {0..}"
    by (intro has_integral_cong) (auto simp: powr_minus field_simps) 
  finally have int:  .  
  have "exp (1 / 2 :: real) ^ 2  2 ^ 2"
    using exp_le by (subst exp_double [symmetric]) simp_all
  hence exp_half: "exp (1 / 2 :: real)  2"
    by (rule power2_le_imp_le) auto

  have mono: "ln x / x ^ 2  ln y / y ^ 2" if "y  exp (1/2)" "x  y" for x y :: real
  proof (rule DERIV_nonpos_imp_nonincreasing[of _ _ "λx. ln x / x ^ 2"])
    fix t assume t: "t  y" "t  x"
    have "y > 0" by (rule less_le_trans[OF _ that(1)]) auto
    with t that have "ln t  ln (exp (1 / 2))"
      by (subst ln_le_cancel_iff) auto
    hence "ln t  1 / 2" by (simp only: ln_exp)
    from t y > 0 have "((λx. ln x / x ^ 2) has_field_derivative ((1 - 2 * ln t) / t ^ 3)) (at t)"
      by (auto intro!: derivative_eq_intros simp: eval_nat_numeral field_simps)
    moreover have "(1 - 2 * ln t) / t ^ 3  0"
      using t that y > 0 ln t  1 / 2 by (intro divide_nonpos_pos) auto
    ultimately show "f'. ((λx. ln x / x ^ 2) has_field_derivative f') (at t)  f'  0" by blast
  qed fact+

  have "fds_converges (fds_deriv fds_zeta) (2 :: complex)"
    by (intro fds_converges_deriv) auto
  hence "(λn. of_real (-ln (real (Suc n)) / (of_nat (Suc n)) ^ 2)) sums deriv zeta 2"
    by (auto simp: fds_converges_altdef add_ac eval_fds_deriv_zeta fds_nth_deriv scaleR_conv_of_real 
             simp del: of_nat_Suc)
  note * = sums_split_initial_segment[OF sums_minus[OF sums_Re[OF this]], of 3]
    have "(λn. ln (real (n+4)) / real (n+4) ^ 2) sums (-Re (deriv zeta 2) - (ln 2 / 4 + ln 3 / 9))"
      using * by (simp add: eval_nat_numeral)
  hence "-Re (deriv zeta 2) - (ln 2 / 4 + ln 3 / 9) =
            (n. ln (real (Suc n) + 3) / (real (Suc n) + 3) ^ 2)"
    by (simp_all add: sums_iff algebra_simps)
  also have "  (ln 3 + 1) / 3" using int exp_half
    by (intro decreasing_sum_le_integral divide_nonneg_pos mono) (auto simp: powr_minus field_simps)
  finally have "-Re (deriv zeta 2)  (16 * ln 3 + 9 * ln 2 + 12) / 36"
    by simp
  also have "ln 3  (11 / 10 :: real)"
    using ln_approx_bounds[of 3 2] by (simp add: power_numeral_reduce numeral_2_eq_2)
  hence "(16 * ln 3 + 9 * ln 2 + 12) / 36  (16 * (11 / 10) + 9 * 25 / 36 + 12) / (36 :: real)"
    using ln2_le_25_over_36 by (intro add_mono mult_left_mono divide_right_mono) auto
  also have " < 1" by simp
  finally show ?thesis by simp
qed

text ‹
  Using the logarithmic derivative of Euler's product formula for $\zeta(s)$ at $s = 2$
  and the bound on $\zeta'(2)$ we have just derived, we can obtain the bound
  \[\sum_{p^i \leq x, i \geq 2} \frac{\ln p}{p^i} < \frac{9}{\pi^2}\ .\]
›
lemma mertens_remainder_aux_bound:
  fixes x :: real
  defines "R  ((p,i) | prime p  i > 1  real (p ^ i)  x. ln (real p) / p ^ i)"
  shows   "R < 9 / pi2"
proof -
  define S' where "S' = {(p, i). prime p  i > 1  real (p ^ i)  x}"
  define S'' where "S'' = {(p, i). prime p  i > 1  real (p ^ Suc i)  x}"

  have finite_row: "finite {i. i > 1  real (p ^ (i + k))  x}" if p: "prime p" for p k
  proof (rule finite_subset)
    show "{i. i > 1  real (p ^ (i + k))  x}  {..nat x}"
    proof safe
      fix i assume i: "i > 1" "real (p ^ (i + k))  x"
      have "i < 2 ^ (i + k)" by (induction i) auto
      also from p have "  p ^ (i + k)" by (intro power_mono) (auto dest: prime_gt_1_nat) 
      also have "real   x" using i by simp
      finally show "i  nat x" by linarith
    qed
  qed auto

  have "S''  S'" unfolding S''_def S'_def
  proof safe
    fix p i assume pi: "prime p" "real (p ^ Suc i)  x" "i > 1"
    have "real (p ^ i)  real (p ^ Suc i)"
      using pi unfolding of_nat_le_iff by (intro power_increasing) (auto dest: prime_gt_1_nat)
    also have "  x" by fact
    finally show "real (p ^ i)  x" .
  qed

  have S'_alt:  "S' = (SIGMA p:{p. prime p  real p  x}. {i. i > 1  real (p ^ i)  x})"
    unfolding S'_def
  proof safe
    fix p i assume "prime p" "real (p ^ i)  x" "i > 1"
    hence "p ^ 1  p ^ i"
      by (intro power_increasing) (auto dest: prime_gt_1_nat)
    also have "real   x" by fact
    finally show "real p  x" by simp
  qed

  have finite: "finite {p. prime p  real p  x}"
    by (rule finite_subset[OF _ finite_Nats_le_real[of x]]) (auto dest: prime_gt_0_nat)
  have "finite S'" unfolding S'_alt using finite_row[of _ 0]
    by (intro finite_SigmaI finite) auto

  have "R  3 / 2 * ((p, i) | (p, i)  S'  even i. ln (real p) / real (p ^ i))"
  proof -
    have "R = (y{0, 1}. z | z  S'  snd z mod 2 = y. ln (real (fst z)) / real (fst z ^ snd z))"
      using finite S' by (subst sum.group) (auto simp: case_prod_unfold R_def S'_def)
    also have " = ((p,i) | (p, i)  S'  even i. ln (real p) / real (p ^ i)) +
                    ((p,i) | (p, i)  S'  odd i. ln (real p) / real (p ^ i))"
      unfolding even_iff_mod_2_eq_zero odd_iff_mod_2_eq_one by (simp add: case_prod_unfold)
    also have "((p,i) | (p, i)  S'  odd i. ln (real p) / real (p ^ i)) =
                 ((p,i) | (p, i)  S''  even i. ln (real p) / real (p ^ Suc i))"
      by (intro sum.reindex_bij_witness[of _ "λ(p,i). (p, Suc i)" "λ(p,i). (p, i - 1)"])
         (auto simp: case_prod_unfold S'_def S''_def elim: oddE simp del: power_Suc)
    also have "  ((p,i) | (p, i)  S'  even i. ln (real p) / real (p ^ Suc i))"
      using S''  S' unfolding case_prod_unfold
      by (intro sum_mono2 divide_nonneg_pos ln_ge_zero finite_subset[OF _ finite S'])
         (auto simp: S'_def S''_def case_prod_unfold dest: prime_gt_0_nat simp del: power_Suc)
    also have "  ((p,i) | (p, i)  S'  even i. ln (real p) / real (2 * p ^ i))"
      unfolding case_prod_unfold
      by (intro sum_mono divide_left_mono) (auto simp: S'_def dest!: prime_gt_1_nat)
    also have " = (1 / 2) * ((p,i) | (p, i)  S'  even i. ln (real p) / real (p ^ i))"
      by (subst sum_distrib_left) (auto simp: case_prod_unfold)
    also have "((p,i) | (p, i)  S'  even i. ln (real p) / real (p ^ i)) +  =
                   3 / 2 * ((p,i) | (p, i)  S'  even i. ln (real p) / real (p ^ i))"
      by simp
    finally show ?thesis by simp
  qed

  also have "((p,i) | (p, i)  S'  even i. ln (real p) / real (p ^ i)) =
               (p | prime p  real p  x. ln (real p) * 
                 (i | i > 0  even i  real (p ^ i)  x. (1 / real p) ^ i))"
   unfolding sum_distrib_left
  proof (subst sum.Sigma[OF _ ballI])
    fix p assume p: "p  {p. prime p  real p  x}"
    thus "finite {i. 0 < i  even i  real (p ^ i)  x}"
      by (intro finite_subset[OF _ exponents_le_finite[of p 1 0 x]]) (auto dest: prime_gt_1_nat)
  qed (auto intro!: sum.cong finite_subset[OF _ finite_Nats_le_real[of x]]
            dest: prime_gt_0_nat simp: S'_alt power_divide)
  also have "  (p | prime p  real p  x. ln (real p) / (real p ^ 2 - 1))"
  proof (rule sum_mono)
    fix p assume p: "p  {p. prime p  real p  x}"
    have "p > 1" using p by (auto dest: prime_gt_1_nat)
    have "(i | i > 0  even i  real (p ^ i)  x. (1 / real p) ^ i) =
            (i | real (p ^ (2 * i + 2))  x. (1 / real p) ^ (2 * i)) / real p ^ 2"
      (is "_ = ?S / _") unfolding sum_divide_distrib
      by (rule sum.reindex_bij_witness[of _ "λi. 2 * Suc i" "λi. (i - 2) div 2"])
         (insert p > 1, auto simp: numeral_3_eq_3 power2_eq_square power_diff
                 algebra_simps elim!: evenE)
    also have "?S = (i | real (p ^ (2 * i + 2))  x. (1 / real p ^ 2) ^ i)"
      by (subst power_mult) (simp_all add: algebra_simps power_divide)
    also have "  (i. (1 / real p ^ 2) ^ i)"
      using exponents_le_finite[of p 2 2 x] p > 1
      by (intro sum_le_suminf) (auto simp: summable_geometric_iff)
    also have " = real p ^ 2 / (real p ^ 2 - 1)"
      using p > 1 by (subst suminf_geometric) (auto simp: field_simps)
    also have " / real p ^ 2 = 1 / (real p ^ 2 - 1)"
      using p > 1 by (simp add: divide_simps)
    finally have "(i | 0 < i  even i  real (p ^ i)  x. (1 / real p) ^ i) 
                     1 / (real p ^ 2 - 1)" (is "?lhs  ?rhs")
      using p > 1 by (simp add: divide_right_mono)
    thus "ln (real p) * ?lhs  ln (real p) / (real p ^ 2 - 1)"
      using p > 1 by (simp add: divide_simps)
  qed
  also have " = (a p | prime p  real p  x. ln (real p) / (real p ^ 2 - 1))"
    using finite by (intro infsetsum_finite [symmetric]) auto
  also have "  (a p | prime p. ln (real p) / (real p ^ 2 - 1))"
    using eval_fds_logderiv_zeta_real[of 2] finite
    by (intro infsetsum_mono_neutral_left divide_nonneg_pos) (auto simp: dest: prime_gt_1_nat)
  also have " = -Re (deriv zeta (of_real 2) / zeta (of_real 2))"
    by (subst eval_fds_logderiv_zeta_real) auto
  also have " = (-Re (deriv zeta 2)) * (6 / pi2)"
    by (simp add: zeta_even_numeral)
  also have " < 1 * (6 / pi2)"
    using deriv_zeta_2_bound by (intro mult_strict_right_mono) auto
  also have "3 / 2 *  = 9 / pi2" by simp
  finally show ?thesis by simp
qed

text ‹
  We now consider the equation
  \[\ln (n!) = \sum_{k\leq n} \Lambda(k) \left\lfloor\frac{n}{k}\right\rfloor\]
  and estimate both sides in different ways. The left-hand-side can be estimated
  using Stirling's formula, and we can simplify the right-hand side to
  \[\sum_{k\leq n} \Lambda(k) \left\lfloor\frac{n}{k}\right\rfloor =
      \sum_{p^i \leq x, i \geq 1} \ln p \left\lfloor\frac{n}{p^i}\right\rfloor\]
  and then split the sum into those $p^i$ with $i = 1$ and those with $i \geq 2$.
  Applying the bound we have just shown and some more routine estimates, we obtain the
  following reasonably strong version of Mertens' First Theorem on the naturals:
  $\mathfrak M(n) - \ln(n) \in (-1-\frac{9}{\pi^2}; \ln 4]$
›
theorem mertens_bound_strong:
  fixes n :: nat assumes n: "n > 0"
  shows   "𝔐 n - ln n  {-1 - 9 / pi2<..ln 4}"
proof (cases "n  3")
  case False
  with n consider "n = 1" | "n = 2" by force
  thus ?thesis
  proof cases
    assume [simp]: "n = 1"
    have "-1 + (-9 / pi2) < 0"
      by (intro add_neg_neg divide_neg_pos) auto
    thus ?thesis by simp
  next
    assume [simp]: "n = 2"
    have eq: "𝔐 n - ln n = -ln 2 / 2" by (simp add: eval_𝔐)
    have "-1 - 9 / pi ^ 2 + ln 2 / 2  -1 - 9 / 4 ^ 2 + 25 / 36 / 2"
      using pi_less_4 ln2_le_25_over_36
      by (intro diff_mono add_mono divide_left_mono divide_right_mono power_mono) auto
    also have " < 0" by simp
    finally have "-ln 2 / 2 > -1 - 9 / pi2" by simp
    moreover {
      have "-ln 2 / 2  (0::real)" by (intro divide_nonpos_pos) auto
      also have "  ln 4" by simp
      finally have "-ln 2 / 2  ln (4 :: real)" by simp
    }
    ultimately show ?thesis unfolding eq by simp
  qed

next
  case True
  hence n: "n  3" by simp
  have finite: "finite {(p, i). prime p  i  1  p ^ i  n}"
  proof (rule finite_subset)
    show "{(p, i). prime p  i  1  p ^ i  n}
             {..nat root 1 (real n)} × {..nat log 2 (real n)}"
      using primepows_le_subset[of "real n" 1] n unfolding of_nat_le_iff by auto
  qed auto

  define r where "r = prime_sum_upto (λp. ln (real p) * frac (real n / real p)) n"
  define R where "R = ((p,i) | prime p  i > 1  p ^ i  n. ln (real p) * real (n div (p ^ i)))"
  define R' where "R' = ((p,i) | prime p  i > 1  p ^ i  n. ln (real p) / p ^ i)"
  have [simp]: "ln (4 :: real) = 2 * ln 2"
    using ln_realpow[of 2 2] by simp
  from pi_less_4 have "ln pi  ln 4" by (subst ln_le_cancel_iff) auto
  also have " = 2 * ln 2" by simp
  also have "  2 * (25 / 36)" by (intro mult_left_mono ln2_le_25_over_36) auto
  finally have ln_pi: "ln pi  25 / 18" by simp
  have "ln 3  ln (4::nat)" by (subst ln_le_cancel_iff) auto
  also have " = 2 * ln 2" by simp
  also have "  2 * (25 / 36)" by (intro mult_left_mono ln2_le_25_over_36) auto
  finally have ln_3: "ln (3::real)  25 / 18" by simp
  
  have "R / n = ((p,i) | prime p  i > 1  p ^ i  n. ln (real p) * (real (n div (p ^ i)) / n))"
    by (simp add: R_def sum_divide_distrib field_simps case_prod_unfold)
  also have "  ((p,i) | prime p  i > 1  p ^ i  n. ln (real p) * (1 / p ^ i))"
    unfolding R'_def case_prod_unfold using n
    by (intro sum_mono mult_left_mono) (auto simp: field_simps real_of_nat_div dest: prime_gt_0_nat)
  also have " = R'" by (simp add: R'_def)
  also have "R' < 9 / pi2"
    unfolding R'_def using mertens_remainder_aux_bound[of n] by simp
  finally have "R / n < 9 / pi2" .
  moreover have "R  0"
    unfolding R_def by (intro sum_nonneg mult_nonneg_nonneg) (auto dest: prime_gt_0_nat)
  ultimately have R_bounds: "R / n  {0..<9 / pi2}" by simp

  have "ln (fact n :: real)  ln (2 * pi * n) / 2 + n * ln n - n + 1 / (12 * n)"
    using ln_fact_bounds(2)[of n] n by simp
  also have " / n - ln n = -1 + (ln 2 + ln pi) / (2 * n) + (ln n / n) / 2 + 1 / (12 * real n ^ 2)"
    using n by (simp add: power2_eq_square field_simps ln_mult)
  also have "  -1 + (ln 2 + ln pi) / (2 * 3) + (ln 3 / 3) / 2 + 1 / (12 * 32)"
    using exp_le n pi_gt3
    by (intro add_mono divide_right_mono divide_left_mono mult_mono
              mult_pos_pos ln_x_over_x_mono power_mono) auto
  also have "  -1 + (25 / 36 + 25 / 18) / (2 * 3) + (25 / 18 / 3) / 2 + 1 / (12 * 32)"
    using ln_pi ln2_le_25_over_36 ln_3 by (intro add_mono divide_left_mono divide_right_mono) auto
  also have "  0" by simp
  finally have "ln n - ln (fact n) / n  0" using n by (simp add: divide_right_mono)
  have "-ln (fact n)  -ln (2 * pi * n) / 2 - n * ln n + n"
    using ln_fact_bounds(1)[of n] n by simp
  also have "ln n +  / n = -ln (2 * pi) / (2 * n) - (ln n / n) / 2 + 1"
    using n by (simp add: field_simps ln_mult)
  also have "  0 - 0 + 1"
    using pi_gt3 n by (intro add_mono diff_mono) auto
  finally have upper: "ln n - ln (fact n) / n  1"
    using n by (simp add: divide_right_mono)
  with ln n - ln (fact n) / n  0 have fact_bounds: "ln n - ln (fact n) / n  {0..1}" by simp

  have "r  prime_sum_upto (λp. ln p * 1) n"
    using less_imp_le[OF frac_lt_1] unfolding r_def θ_def prime_sum_upto_def
    by (intro sum_mono mult_left_mono) (auto simp: dest: prime_gt_0_nat)
  also have " = θ n" by (simp add: θ_def)
  also have " < ln 4 * n" using n by (intro θ_upper_bound) auto
  finally have "r / n < ln 4" using n by (simp add: field_simps)
  moreover have "r  0" unfolding r_def prime_sum_upto_def
    by (intro sum_nonneg mult_nonneg_nonneg) (auto dest: prime_gt_0_nat)
  ultimately have r_bounds: "r / n  {0..<ln 4}" by simp

  have "ln (fact n :: real) = sum_upto (λk. mangoldt k * real (n div k)) (real n)"
    by (simp add: ln_fact_conv_sum_upto_mangoldt)
  also have " = ((p,i) | prime p  i > 0  real (p ^ i)  real n.
                       ln (real p) * real (n div (p ^ i)))"
    by (intro sum_upto_primepows) (auto simp: mangoldt_non_primepow)
  also have "{(p, i). prime p  i > 0  real (p ^ i)  real n} =
               {(p, i). prime p  i = 1  p  n} 
               {(p, i). prime p  i > 1  (p ^ i)  n}" unfolding of_nat_le_iff
    by (auto simp: not_less le_Suc_eq)
  also have "((p,i). ln (real p) * real (n div (p ^ i))) =
               ((p,i) | prime p  i = 1  p  n. ln (real p) * real (n div (p ^ i))) + R"
    (is "_ = ?S + _")
    by (subst sum.union_disjoint) (auto intro!: finite_subset[OF _ finite] simp: R_def)
  also have "?S = prime_sum_upto (λp. ln (real p) * real (n div p)) n"
    unfolding prime_sum_upto_def
    by (intro sum.reindex_bij_witness[of _ "λp. (p, 1)" fst]) auto
  also have " = prime_sum_upto (λp. ln (real p) * real n / real p) n - r"
    unfolding r_def prime_sum_upto_def sum_subtractf[symmetric] using n
    by (intro sum.cong) (auto simp: frac_def real_of_nat_div algebra_simps)
  also have "prime_sum_upto (λp. ln (real p) * real n / real p) n = n * 𝔐 n"
    by (simp add: primes_M_def sum_distrib_left sum_distrib_right prime_sum_upto_def field_simps)
  finally have "𝔐 n - ln n = ln (fact n) / n - ln n + r / n - R / n"
    using n by (simp add: field_simps)
  hence "ln n - 𝔐 n = ln n - ln (fact n) / n - r / n + R / n"
    by simp
  with fact_bounds r_bounds R_bounds show "𝔐 n - ln n  {-1 - 9 / pi2<..ln 4}"
    by simp
qed

text ‹
  As a simple corollary, we obtain a similar bound on the reals.
›
lemma mertens_bound_real_strong:
  fixes x :: real assumes x: "x  1"
  shows   "𝔐 x - ln x  {-1 - 9 / pi ^ 2 - ln (1 + frac x / real (nat x)) <.. ln 4}"
proof -
  have "𝔐 x - ln x  𝔐 (real (nat x)) - ln (real (nat x))"
    using assms by simp
  also have "  ln 4"
    using mertens_bound_strong[of "nat x"] assms by simp
  finally have "𝔐 x - ln x  ln 4" .

  from assms have pos: "real_of_int x  0" by linarith
  have "frac x / real (nat x)  0"
    using assms by (intro divide_nonneg_pos) auto
  moreover have "frac x / real (nat x)  1 / 1"
    using assms frac_lt_1[of x] by (intro frac_le) auto
  ultimately have *: "frac x / real (nat x)  {0..1}" by auto
  have "ln x - ln (real (nat x)) = ln (x / real (nat x))"
    using assms by (subst ln_div) auto
  also have "x / real (nat x) = 1 + frac x / real (nat x)"
    using assms pos by (simp add: frac_def field_simps)
  finally have "𝔐 x - ln x > -1-9/pi^2-ln (1 + frac x / real (nat x))"
    using mertens_bound_strong[of "nat x"] x by simp
  with 𝔐 x - ln x  ln 4 show ?thesis by simp
qed

text ‹
  We weaken this estimate a bit to obtain nicer bounds:
›
lemma mertens_bound_real':
  fixes x :: real assumes x: "x  1"
  shows   "𝔐 x - ln x  {-2<..25/18}"
proof -
  have "𝔐 x - ln x  ln 4"
    using mertens_bound_real_strong[of x] x by simp
  also have "  25 / 18"
    using ln_realpow[of 2 2] ln2_le_25_over_36 by simp
  finally have "𝔐 x - ln x  25 / 18" .

  have ln2: "ln (2 :: real)  {2/3..25/36}"
    using ln_approx_bounds[of 2 1] by (simp add: eval_nat_numeral)
  have ln3: "ln (3::real)  {1..10/9}"
    using ln_approx_bounds[of 3 1] by (simp add: eval_nat_numeral)
  have ln5: "ln (5::real)  {4/3..76/45}"
    using ln_approx_bounds[of 5 1] by (simp add: eval_nat_numeral)
  have ln7: "ln (7::real)  {3/2..15/7}"
    using ln_approx_bounds[of 7 1] by (simp add: eval_nat_numeral)
  have ln11: "ln (11::real)  {5/3..290/99}"
    using ln_approx_bounds[of 11 1] by (simp add: eval_nat_numeral)

  ― ‹Choosing the lower bound -2 is somewhat arbitrary here; it is a trade-off between getting
      a reasonably tight bound and having to make lots of case distinctions. To get -2 as a lower
      bound, we have to show the cases up to x = 11› by case distinction,›
  have "𝔐 x - ln x > -2"
  proof (cases "x  11")
    case False
    hence "x  {1..<2}  x  {2..<3}  x  {3..<5}  x  {5..<7}  x  {7..<11}"
      using x by force
    thus ?thesis
    proof (elim disjE)
      assume x: "x  {1..<2}"
      hence "ln x - 𝔐 x  ln 2 - 0"
        by (intro diff_mono) auto
      also have " < 2" using ln2_le_25_over_36 by simp
      finally show ?thesis by simp
    next
      assume x: "x  {2..<3}"
      hence [simp]: "x = 2" by (intro floor_unique) auto
      from x have "ln x - 𝔐 x  ln 3 - ln 2 / 2"
        by (intro diff_mono) (auto simp: eval_𝔐)
      also have " = ln (9 / 2) / 2" using ln_realpow[of 3 2] by (simp add: ln_div)
      also have " < 2" using ln_approx_bounds[of "9 / 2" 1] by (simp add: eval_nat_numeral)
      finally show ?thesis by simp
    next
      assume x: "x  {3..<5}"
      hence "𝔐 3 = 𝔐 x"
        unfolding primes_M_def
        by (intro prime_sum_upto_eqI'[where a' = 3 and b' = 4])
           (auto simp: nat_le_iff le_numeral_iff nat_eq_iff floor_eq_iff)
      also have "𝔐 3 = ln 2 / 2 + ln 3 / 3"
        by (simp add: eval_𝔐 eval_nat_numeral mark_out_code)
      finally have [simp]: "𝔐 x = ln 2 / 2 + ln 3 / 3" ..
      from x have "ln x - 𝔐 x  ln 5 - (ln 2 / 2 + ln 3 / 3)"
        by (intro diff_mono) auto
      also have " < 2" using ln2 ln3 ln5 by simp
      finally show ?thesis by simp
    next
      assume x: "x  {5..<7}"
      hence "𝔐 5 = 𝔐 x"
        unfolding primes_M_def
        by (intro prime_sum_upto_eqI'[where a' = 5 and b' = 6])
           (auto simp: nat_le_iff le_numeral_iff nat_eq_iff floor_eq_iff)
      also have "𝔐 5 = ln 2 / 2 + ln 3 / 3 + ln 5 / 5"
        by (simp add: eval_𝔐 eval_nat_numeral mark_out_code)
      finally have [simp]: "𝔐 x = ln 2 / 2 + ln 3 / 3 + ln 5 / 5" ..
      from x have "ln x - 𝔐 x  ln 7 - (ln 2 / 2 + ln 3 / 3 + ln 5 / 5)"
        by (intro diff_mono) auto
      also have " < 2" using ln2 ln3 ln5 ln7 by simp
      finally show ?thesis by simp
    next
      assume x: "x  {7..<11}"
      hence "𝔐 7 = 𝔐 x"
        unfolding primes_M_def
        by (intro prime_sum_upto_eqI'[where a' = 7 and b' = 10])
           (auto simp: nat_le_iff le_numeral_iff nat_eq_iff floor_eq_iff)
      also have "𝔐 7 = ln 2 / 2 + ln 3 / 3 + ln 5 / 5 + ln 7 / 7"
        by (simp add: eval_𝔐 eval_nat_numeral mark_out_code)
      finally have [simp]: "𝔐 x = ln 2 / 2 + ln 3 / 3 + ln 5 / 5 + ln 7 / 7" ..
      from x have "ln x - 𝔐 x  ln 11 - (ln 2 / 2 + ln 3 / 3 + ln 5 / 5 + ln 7 / 7)"
        by (intro diff_mono) auto
      also have " < 2" using ln2 ln3 ln5 ln7 ln11 by simp
      finally show ?thesis by simp
    qed
  next
    case True
    have "ln x - 𝔐 x  1 + 9/pi^2 + ln (1 + frac x / real (nat x))"
      using mertens_bound_real_strong[of x] x by simp
    also have "1 + frac x / real (nat x)  1 + 1 / 11"
      using True frac_lt_1[of x] by (intro add_mono frac_le) auto
    hence "ln (1 + frac x / real (nat x))  ln (1 + 1 / 11)"
      using x by (subst ln_le_cancel_iff) (auto intro!: add_pos_nonneg)
    also have " = ln (12 / 11)" by simp
    also have "  1585 / 18216"
      using ln_approx_bounds[of "12 / 11" 1] by (simp add: eval_nat_numeral)
    also have "9 / pi ^ 2  9 / 3.141592653588 ^ 2"
      using pi_approx by (intro divide_left_mono power_mono mult_pos_pos) auto
    also have "1 +  + 1585 / 18216 < 2"
      by (simp add: power2_eq_square)
    finally show ?thesis by simp
  qed
  with 𝔐 x - ln x  25 / 18 show ?thesis by simp
qed

corollary mertens_first_theorem:
  fixes x :: real assumes x: "x  1"
  shows   "¦𝔐 x - ln x¦ < 2"
  using mertens_bound_real'[of x] x by (simp add: abs_if)


subsection ‹Mertens' Second Theorem›

text ‹
  Mertens' Second Theorem concerns the asymptotics of the Prime Harmonic Series, namely
  \[\sum\limits_{p \leq x} \frac{1}{p} = \ln\ln x + M + O\left(\frac{1}{\ln x}\right)\]
  where $M \approx 0.261497$ is the Meissel--Mertens constant.

  We define the constant in the following way:
›
definition meissel_mertens where
  "meissel_mertens = 1 - ln (ln 2) + integral {2..} (λt. (𝔐 t - ln t) / (t * ln t ^ 2))"

text ‹
  We will require the value of the integral
  $\int_a^\infty \frac{t}{\ln^2 t} \textrm{d}t = \frac{1}{\ln a}$
  as an upper bound on the remainder term:
›
lemma integral_one_over_x_ln_x_squared:
  assumes a: "(a::real) > 1"
  shows "set_integrable lborel {a<..} (λt. 1 / (t * ln t ^ 2))" (is ?th1)
    and "set_lebesgue_integral lborel {a<..} (λt. 1 / (t * ln t ^ 2)) = 1 / ln a" (is ?th2)
    and "((λt. 1 / (t * (ln t)2)) has_integral 1 / ln a) {a<..}" (is ?th3)
proof -
  have cont: "isCont (λt. 1 / (t * (ln t)2)) x" if "x > a" for x
    using that a by (auto intro!: continuous_intros)
  have deriv: "((λx. -1 / ln x) has_real_derivative 1 / (x * (ln x)2)) (at x)" if "x > a" for x
    using that a by (auto intro!: derivative_eq_intros simp: power2_eq_square field_simps)
  have lim1: "(((λx. - 1 / ln x)  real_of_ereal)  -(1 / ln a)) (at_right (ereal a))"
    unfolding ereal_tendsto_simps using a by (real_asymp simp: field_simps)
  have lim2: "(((λx. - 1 / ln x)  real_of_ereal)  0) (at_left )"
    unfolding ereal_tendsto_simps using a by (real_asymp simp: field_simps)

  have "set_integrable lborel (einterval a ) (λt. 1 / (t * ln t ^ 2))"
    by (rule interval_integral_FTC_nonneg[OF _ deriv cont _ lim1 lim2]) (use a in auto)
  thus ?th1 by simp
  have "interval_lebesgue_integral lborel (ereal a)  (λt. 1 / (t * ln t ^ 2)) = 0 - (-(1 / ln a))"
    by (rule interval_integral_FTC_nonneg[OF _ deriv cont _ lim1 lim2]) (use a in auto)
  thus ?th2 by (simp add: interval_integral_to_infinity_eq)

  have "((λt. 1 / (t * ln t ^ 2)) has_integral
           set_lebesgue_integral lebesgue {a<..} (λt. 1 / (t * ln t ^ 2))) {a<..}"
    using ?th1 by (intro has_integral_set_lebesgue)
                    (auto simp: set_integrable_def integrable_completion)
  also have "set_lebesgue_integral lebesgue {a<..} (λt. 1 / (t * ln t ^ 2)) = 1 / ln a"
    using ?th2 unfolding set_lebesgue_integral_def by (subst integral_completion) auto
  finally show ?th3 .
qed

text ‹
  We show that the integral in our definition of the Meissel--Mertens constant is
  well-defined and give an upper bound for its tails:
›
lemma
  assumes "a > (1 :: real)"
  defines "r  (λt. (𝔐 t - ln t) / (t * ln t ^ 2))"
  shows integrable_meissel_mertens:  "set_integrable lborel {a<..} r"
    and meissel_mertens_integral_le: "norm (integral {a<..} r)  2 / ln a"
proof -
  have *: "((λt. 2 * (1 / (t * ln t ^ 2))) has_integral 2 * (1 / ln a)) {a<..}"
    using assms by (intro has_integral_mult_right integral_one_over_x_ln_x_squared) auto
  show "set_integrable lborel {a<..} r" unfolding set_integrable_def
  proof (rule Bochner_Integration.integrable_bound[OF _ _ AE_I2])
    have "integrable lborel (λt::real. indicator {a<..} t * (2 * (1 / (t * ln t ^ 2))))"
      using integrable_mult_right[of 2, 
              OF integral_one_over_x_ln_x_squared(1)[of a, unfolded set_integrable_def]] assms
      by (simp add: algebra_simps)
    thus "integrable lborel (λt::real. indicator {a<..} t *R (2 / (t * ln t ^ 2)))"
      by simp
    fix x :: real
    show "norm (indicat_real {a<..} x *R r x) 
            norm (indicat_real {a<..} x *R (2 / (x * ln x ^ 2)))"
    proof (cases "x > a")
      case True
      thus ?thesis
      unfolding norm_scaleR norm_mult r_def norm_divide using mertens_first_theorem[of x] assms
        by (intro mult_mono frac_le divide_nonneg_pos) (auto simp: indicator_def)
    qed (auto simp: indicator_def)
  qed (auto simp: r_def)
  hence "r integrable_on {a<..}"
    by (simp add: set_borel_integral_eq_integral(1))
  hence "norm (integral {a<..} r)  integral {a<..} (λx. 2 * (1 / (x * ln x ^ 2)))"
  proof (rule integral_norm_bound_integral)
    show  "(λx. 2 * (1 / (x * (ln x)2))) integrable_on {a<..}"
      using * by (simp add: has_integral_iff)
    fix x assume "x  {a<..}"
    hence "norm (r x)  2 / (x * (ln x)2)"
      unfolding r_def norm_divide using mertens_first_theorem[of x] assms
      by (intro mult_mono frac_le divide_nonneg_pos) (auto simp: indicator_def)
    thus "norm (r x)  2* (1 / (x * ln x ^ 2))" by simp
  qed
  also have " = 2 / ln a"
    using * by (simp add: has_integral_iff)
  finally show "norm (integral {a<..} r)  2 / ln a" .
qed

lemma integrable_on_meissel_mertens:
  assumes "A  {1..}" "Inf A > 1" "A  sets borel"
  shows   "(λt. (𝔐 t - ln t) / (t * ln t ^ 2)) integrable_on A"
proof -
  from assms obtain x where x: "1 < x" "x < Inf A"
    using dense by blast
  from assms have "bdd_below A" by (intro bdd_belowI[of _ 1]) auto
  hence "A  {Inf A..}" by (auto simp: cInf_lower)
  also have "  {x<..}" using x by auto
  finally have *: "A  {x<..}" .
  have "set_integrable lborel A (λt. (𝔐 t - ln t) / (t * ln t ^ 2))"
    by (rule set_integrable_subset[OF integrable_meissel_mertens[of x]]) (use x * assms in auto)
  thus ?thesis by (simp add: set_borel_integral_eq_integral(1))
qed

lemma meissel_mertens_bounds: "¦meissel_mertens - 1 + ln (ln 2)¦  2 / ln 2"
proof -
  have *: "{2..} - {2<..} = {2::real}" by auto
  also have "negligible " by simp
  finally have "integral {2..} (λt. (𝔐 t - ln t) / (t * (ln t)2)) =
                  integral {2<..} (λt. (𝔐 t - ln t) / (t * (ln t)2))"
    by (intro sym[OF integral_subset_negligible]) auto
  also have "norm   2 / ln 2"
    by (rule meissel_mertens_integral_le) auto
  finally show "¦meissel_mertens - 1 + ln (ln 2)¦  2 / ln 2"
    by (simp add: meissel_mertens_def)
qed

text ‹
  Finally, obtaining Mertens' second theorem from the first one is nothing but a routine
  summation by parts, followed by a use of the above bound:
›
theorem mertens_second_theorem:
  defines "f  prime_sum_upto (λp. 1 / p)"
  shows   "x. x  2  ¦f x - ln (ln x) - meissel_mertens¦  4 / ln x"
    and   "(λx. f x - ln (ln x) - meissel_mertens)  O(λx. 1 / ln x)"
proof -
  define r where "r = (λt. (𝔐 t - ln t) / (t * ln t ^ 2))"

  {
    fix x :: real assume x: "x > 2"
    have "((λt. 𝔐 t * (-1 / (t * ln t ^ 2))) has_integral 𝔐 x * (1 / ln x) - 𝔐 2 * (1 / ln 2) -
            (nreal -` {2<..x}. ind prime n * (ln n / real n) * (1 / ln n))) {2..x}"
      unfolding primes_M_def prime_sum_upto_altdef1 using x
      by (intro partial_summation_strong[of "{}"])
         (auto intro!: continuous_intros derivative_eq_intros simp: power2_eq_square
               simp flip: has_real_derivative_iff_has_vector_derivative)
    also have "𝔐 x * (1 / ln x) - 𝔐 2 * (1 / ln 2) - 
                 (nreal -` {2<..x}. ind prime n * (ln n / n) * (1 / ln n)) =
               𝔐 x / ln x - (ninsert 2 (real -` {2<..x}). ind prime n * (ln n / n) * (1 / ln n))" 
      (is "_ = _ - ?S")
      by (subst sum.insert)
         (auto simp: primes_M_def finite_vimage_real_of_nat_greaterThanAtMost eval_prime_sum_upto)
    also have "?S = f x"
      unfolding f_def prime_sum_upto_altdef1 sum_upto_def using x
      by (intro sum.mono_neutral_cong_left) (auto simp: not_less numeral_2_eq_2 le_Suc_eq)
    finally have "((λt. -𝔐 t / (t * ln t ^ 2)) has_integral (𝔐 x / ln x - f x)) {2..x}"
      by simp
    from has_integral_neg[OF this]
      have "((λt. 𝔐 t / (t * ln t ^ 2)) has_integral (f x - 𝔐 x / ln x)) {2..x}" by simp
    hence "((λt. 𝔐 t / (t * ln t ^ 2) - 1 / (t * ln t)) has_integral
             (f x - 𝔐 x / ln x - (ln (ln x) - ln (ln 2)))) {2..x}" using x
      by (intro has_integral_diff fundamental_theorem_of_calculus)
         (auto simp flip: has_real_derivative_iff_has_vector_derivative
               intro!: derivative_eq_intros)
    also have "?this  (r has_integral (f x - 𝔐 x / ln x - (ln (ln x) - ln (ln 2)))) {2..x}"
      by (intro has_integral_cong) (auto simp: r_def field_simps power2_eq_square)
    finally have  .
  } note integral = this

  define R𝔐 where "R𝔐 = (λx. 𝔐 x - ln x)"
  have 𝔐: "𝔐 x = ln x + R𝔐 x" for x by (simp add: R𝔐_def)
  define I where "I = (λx. integral {x..} r)"
  define C where "C = (1 - ln (ln 2) + I 2)"
  have C_altdef: "C = meissel_mertens"
    by (simp add: I_def r_def C_def meissel_mertens_def)

  show bound: " ¦f x - ln (ln x) - meissel_mertens¦  4 / ln x" if x: "x  2" for x
  proof (cases "x = 2")
    case True
    hence "¦f x - ln (ln x) - meissel_mertens¦ = ¦1 / 2 - ln (ln 2) - meissel_mertens¦"
      by (simp add: f_def eval_prime_sum_upto )
    also have "  2 / ln 2 + 1 / 2"
      using meissel_mertens_bounds by linarith
    also have "  2 / ln 2 + 2 / ln 2" using ln2_le_25_over_36
      by (intro add_mono divide_left_mono) auto
    finally show ?thesis using True by simp
  next
    case False
    hence x: "x > 2" using x by simp
    have "integral {2..x} r + I x = integral ({2..x}  {x..}) r" unfolding I_def r_def using x
      by (intro integral_Un [symmetric] integrable_on_meissel_mertens) (auto simp: max_def r_def)
    also have "{2..x}  {x..} = {2..}" using x by auto
    finally have *: "integral {2..x} r = I 2 - I x" unfolding I_def by simp
    have eq: "f x - ln (ln x) - C = R𝔐 x / ln x - I x"
      using integral[OF x] x by (auto simp: C_def field_simps 𝔐 has_integral_iff *)
    also have "¦¦  ¦R𝔐 x / ln x¦ + norm (I x)"
      unfolding real_norm_def by (rule abs_triangle_ineq4)
    also have "¦R𝔐 x / ln x¦  2 / ¦ln x¦"
      unfolding R𝔐_def abs_divide using mertens_first_theorem[of x] x
      by (intro divide_right_mono) auto
    also have "{x..} - {x<..} = {x}" and "{x<..}  {x..}" by auto
    hence "I x = integral {x<..} r" unfolding I_def
      by (intro integral_subset_negligible [symmetric]) simp_all
    also have "norm   2 / ln x"
      using meissel_mertens_integral_le[of x] x by (simp add: r_def)
    finally show "¦f x - ln (ln x) - meissel_mertens¦  4 / ln x"
      using x by (simp add: C_altdef)
  qed

  have "(λx. f x - ln (ln x) - C)  O(λx. 1 / ln x)"
  proof (intro landau_o.bigI[of 4] eventually_mono[OF eventually_ge_at_top[of 2]])
    fix x :: real assume x: "x  2"
    with bound[OF x] show "norm (f x - ln (ln x) - C)  4 * norm (1 / ln x)"
      by (simp add: C_altdef)
  qed (auto intro!: add_pos_nonneg)
  thus "(λx. f x - ln (ln x) - meissel_mertens)  O(λx. 1 / ln x)"
    by (simp add: C_altdef)
qed

corollary prime_harmonic_asymp_equiv: "prime_sum_upto (λp. 1 / p) ∼[at_top] (λx. ln (ln x))"
proof -
  define f where "f = prime_sum_upto (λp. 1 / p)"
  have "(λx. f x - ln (ln x) - meissel_mertens + meissel_mertens)  o(λx. ln (ln x))"
    unfolding f_def
    by (rule sum_in_smallo[OF landau_o.big_small_trans[OF mertens_second_theorem(2)]]) real_asymp+
  hence "(λx. f x - ln (ln x))  o(λx. ln (ln x))"
    by simp
  thus ?thesis unfolding f_def
    by (rule smallo_imp_asymp_equiv)
qed

text ‹
  As a corollary, we get the divergence of the prime harmonic series.
›
corollary prime_harmonic_diverges: "filterlim (prime_sum_upto (λp. 1 / p)) at_top at_top"
  using asymp_equiv_symI[OF prime_harmonic_asymp_equiv]
  by (rule asymp_equiv_at_top_transfer) real_asymp

(*<*)
unbundle no_prime_counting_notation
(*>*)

end