# Theory 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"]
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)"
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 / pi⇧2"
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 / pi⇧2)"
also have "… < 1 * (6 / pi⇧2)"
using deriv_zeta_2_bound by (intro mult_strict_right_mono) auto
also have "3 / 2 * … = 9 / pi⇧2" 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 / pi⇧2<..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 / pi⇧2) < 0"
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 / pi⇧2" 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 / pi⇧2"
unfolding R'_def using mertens_remainder_aux_bound[of n] by simp
finally have "R / n < 9 / pi⇧2" .
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 / pi⇧2}" 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 * 3⇧2)"
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 * 3⇧2)"
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)"
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 / pi⇧2<..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"
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
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<..}"
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"
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) -
(∑n∈real - {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) -
(∑n∈real - {2<..x}. ind prime n * (ln n / n) * (1 / ln n)) =
𝔐 x / ln x - (∑n∈insert 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
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)"
thus "(λx. f x - ln (ln x) - meissel_mertens) ∈ O(λx. 1 / ln x)"
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`