# Theory Bertrand

(*
File:     Bertrand.thy
Authors:  Julian Biendarra, Manuel Eberl <eberlm@in.tum.de>, Larry Paulson

A proof of Bertrand's postulate (based on John Harrison's HOL Light proof).
Uses reflection and the approximation tactic.
*)
theory Bertrand
imports
Complex_Main
"HOL-Number_Theory.Number_Theory"
"HOL-Library.Discrete"
"HOL-Decision_Procs.Approximation_Bounds"
"HOL-Library.Code_Target_Numeral"
Pratt_Certificate.Pratt_Certificate
begin

subsection ‹Auxiliary facts›

lemma ln_2_le: "ln 2 ≤ 355 / (512 :: real)"
proof -
have "ln 2 ≤ real_of_float (ub_ln2 12)" by (rule ub_ln2)
also have "ub_ln2 12 = Float 5680 (- 13)" by code_simp
finally show ?thesis by simp
qed

lemma ln_2_ge: "ln 2 ≥ (5677 / 8192 :: real)"
proof -
have "ln 2 ≥ real_of_float (lb_ln2 12)" by (rule lb_ln2)
also have "lb_ln2 12 = Float 5677 (-13)" by code_simp
finally show ?thesis by simp
qed

lemma ln_2_ge': "ln (2 :: real) ≥ 2/3" and ln_2_le': "ln (2 :: real) ≤ 16/23"
using ln_2_le ln_2_ge by simp_all

lemma of_nat_ge_1_iff: "(of_nat x :: 'a :: linordered_semidom) ≥ 1 ⟷ x ≥ 1"
using of_nat_le_iff[of 1 x] by (subst (asm) of_nat_1)

lemma floor_conv_div_nat:
"of_int (floor (real m / real n)) = real (m div n)"
by (subst floor_divide_of_nat_eq) simp

lemma frac_conv_mod_nat:
"frac (real m / real n) = real (m mod n) / real n"
by (cases "n = 0")
(simp_all add: frac_def floor_conv_div_nat field_simps of_nat_mult

lemma of_nat_prod_mset: "prod_mset (image_mset of_nat A) = of_nat (prod_mset A)"
by (induction A) simp_all

lemma prod_mset_pos: "(⋀x :: 'a :: linordered_semidom. x ∈# A ⟹ x > 0) ⟹ prod_mset A > 0"
by (induction A) simp_all

lemma ln_msetprod:
assumes "⋀x. x ∈#I ⟹ x > 0"
shows "(∑p::nat∈#I. ln p) = ln (∏p∈#I. p)"
using assms by (induction I) (simp_all add: of_nat_prod_mset ln_mult prod_mset_pos)

lemma ln_fact: "ln (fact n) = (∑d=1..n. ln d)"
by (induction n) (simp_all add: ln_mult)

lemma overpower_lemma:
fixes f g :: "real ⇒ real"
assumes "f a ≤ g a"
assumes "⋀x. a ≤ x ⟹ ((λx. g x - f x) has_real_derivative (d x)) (at x)"
assumes "⋀x. a ≤ x ⟹ d x ≥ 0"
assumes "a ≤ x"
shows   "f x ≤ g x"
proof (cases "a < x")
case True
with assms have "∃z. z > a ∧ z < x ∧ g x - f x - (g a - f a) = (x - a) * d z"
by (intro MVT2) auto
then obtain z where z: "z > a" "z < x" "g x - f x - (g a - f a) = (x - a) * d z" by blast
hence "f x = g x + (f a - g a) + (a - x) * d z" by (simp add: algebra_simps)
also from assms have "f a - g a ≤ 0" by (simp add: algebra_simps)
also from assms z have "(a - x) * d z ≤ 0 * d z"
by (intro mult_right_mono) simp_all
finally show ?thesis by simp
qed (insert assms, auto)

subsection ‹Preliminary definitions›

definition primepow_even :: "nat ⇒ bool" where
"primepow_even q ⟷ (∃ p k. 1 ≤ k ∧ prime p ∧ q = p^(2*k))"

definition primepow_odd :: "nat ⇒ bool" where
"primepow_odd q ⟷ (∃ p k. 1 ≤ k ∧ prime p ∧ q = p^(2*k+1))"

abbreviation (input) isprimedivisor :: "nat ⇒ nat ⇒ bool" where
"isprimedivisor q p ≡ prime p ∧ p dvd q"

definition pre_mangoldt :: "nat ⇒ nat" where
"pre_mangoldt d = (if primepow d then aprimedivisor d else 1)"

definition mangoldt_even :: "nat ⇒ real" where
"mangoldt_even d = (if primepow_even d then ln (real (aprimedivisor d)) else 0)"

definition mangoldt_odd :: "nat ⇒ real" where
"mangoldt_odd d = (if primepow_odd d then ln (real (aprimedivisor d)) else 0)"

definition mangoldt_1 :: "nat ⇒ real" where
"mangoldt_1 d = (if prime d then ln d else 0)"

definition psi :: "nat ⇒ real" where
"psi n = (∑d=1..n. mangoldt d)"

definition psi_even :: "nat ⇒ real" where
"psi_even n = (∑d=1..n. mangoldt_even d)"

definition psi_odd :: "nat ⇒ real" where
"psi_odd n = (∑d=1..n. mangoldt_odd d)"

abbreviation (input) psi_even_2 :: "nat ⇒ real" where
"psi_even_2 n ≡ (∑d=2..n. mangoldt_even d)"

abbreviation (input) psi_odd_2 :: "nat ⇒ real" where
"psi_odd_2 n ≡ (∑d=2..n. mangoldt_odd d)"

definition theta :: "nat ⇒ real" where
"theta n = (∑p=1..n. if prime p then ln (real p) else 0)"

subsection ‹Properties of prime powers›

lemma primepow_even_imp_primepow:
assumes "primepow_even n"
shows   "primepow n"
proof -
from assms obtain p k where "1 ≤ k" "prime p" "n = p ^ (2 * k)"
unfolding primepow_even_def by blast
moreover from ‹1 ≤ k› have "2 * k > 0"
by simp
ultimately show ?thesis unfolding primepow_def by blast
qed

lemma primepow_odd_imp_primepow:
assumes "primepow_odd n"
shows   "primepow n"
proof -
from assms obtain p k where "1 ≤ k" "prime p" "n = p ^ (2 * k + 1)"
unfolding primepow_odd_def by blast
moreover from ‹1 ≤ k› have "Suc (2 * k) > 0"
by simp
ultimately show ?thesis unfolding primepow_def
by (auto simp del: power_Suc)
qed

lemma primepow_odd_altdef:
"primepow_odd n ⟷
primepow n ∧ odd (multiplicity (aprimedivisor n) n) ∧ multiplicity (aprimedivisor n) n > 1"
proof (intro iffI conjI; (elim conjE)?)
assume "primepow_odd n"
then obtain p k where n: "k ≥ 1" "prime p" "n = p ^ (2 * k + 1)"
by (auto simp: primepow_odd_def)
thus "odd (multiplicity (aprimedivisor n) n)" "multiplicity (aprimedivisor n) n > 1"
next
assume A: "primepow n" and B: "odd (multiplicity (aprimedivisor n) n)"
and C: "multiplicity (aprimedivisor n) n > 1"
from A obtain p k where n: "k ≥ 1" "prime p" "n = p ^ k"
by (auto simp: primepow_def Suc_le_eq)
with B C have "odd k" "k > 1"
then obtain j where j: "k = 2 * j + 1" "j > 0" by (auto elim!: oddE)
with n show "primepow_odd n" by (auto simp: primepow_odd_def intro!: exI[of _ p, OF exI[of _ j]])
qed (auto dest: primepow_odd_imp_primepow)

lemma primepow_even_altdef:
"primepow_even n ⟷ primepow n ∧ even (multiplicity (aprimedivisor n) n)"
proof (intro iffI conjI; (elim conjE)?)
assume "primepow_even n"
then obtain p k where n: "k ≥ 1" "prime p" "n = p ^ (2 * k)"
by (auto simp: primepow_even_def)
thus "even (multiplicity (aprimedivisor n) n)"
next
assume A: "primepow n" and B: "even (multiplicity (aprimedivisor n) n)"
from A obtain p k where n: "k ≥ 1" "prime p" "n = p ^ k"
by (auto simp: primepow_def Suc_le_eq)
with B have "even k"
then obtain j where j: "k = 2 * j" by (auto elim!: evenE)
from j n have "j ≠ 0" by (intro notI) simp_all
with j n show "primepow_even n"
by (auto simp: primepow_even_def intro!: exI[of _ p, OF exI[of _ j]])
qed (auto dest: primepow_even_imp_primepow)

lemma primepow_odd_mult:
assumes "d > Suc 0"
shows   "primepow_odd (aprimedivisor d * d) ⟷ primepow_even d"
using assms
by (auto simp: primepow_odd_altdef primepow_even_altdef primepow_mult_aprimedivisorI
aprimedivisor_primepow prime_aprimedivisor' aprimedivisor_dvd'
prime_elem_multiplicity_mult_distrib prime_elem_aprimedivisor_nat
dest!: primepow_multD)

lemma pre_mangoldt_primepow:
assumes "primepow n" "aprimedivisor n = p"
shows   "pre_mangoldt n = p"
using assms by (simp add: pre_mangoldt_def)

lemma pre_mangoldt_notprimepow:
assumes "¬primepow n"
shows   "pre_mangoldt n = 1"
using assms by (simp add: pre_mangoldt_def)

lemma primepow_cases:
"primepow d ⟷
(  primepow_even d ∧ ¬ primepow_odd d ∧ ¬ prime d) ∨
(¬ primepow_even d ∧   primepow_odd d ∧ ¬ prime d) ∨
(¬ primepow_even d ∧ ¬ primepow_odd d ∧   prime d)"
by (auto simp: primepow_even_altdef primepow_odd_altdef multiplicity_aprimedivisor_Suc_0_iff
elim!: oddE intro!: Nat.gr0I)

subsection ‹Deriving a recurrence for the psi function›

lemma ln_fact_bounds:
assumes "n > 0"
shows "abs(ln (fact n) - n * ln n + n) ≤ 1 + ln n"
proof -
have "∀n∈{0<..}. ∃z>real n. z < real (n + 1) ∧ real (n + 1) * ln (real (n + 1)) -
real n * ln (real n) = (real (n + 1) - real n) * (ln z + 1)"
by (intro ballI MVT2) (auto intro!: derivative_eq_intros)
hence "∀n∈{0<..}. ∃z>real n. z < real (n + 1) ∧ real (n + 1) * ln (real (n + 1)) -
real n * ln (real n) = (ln z + 1)" by (simp add: algebra_simps)
from bchoice[OF this] obtain k :: "nat ⇒ real"
where lb: "real n < k n" and ub: "k n < real (n + 1)" and
mvt: "real (n+1) * ln (real (n+1)) - real n * ln (real n) = ln (k n) + 1"
if "n > 0" for n::nat by blast
have *: "(n + 1) * ln (n + 1) = (∑i=1..n. ln(k i) + 1)" for n::nat
proof (induction n)
case (Suc n)
have "(∑i = 1..n+1. ln (k i) + 1) = (∑i = 1..n. ln (k i) + 1) + ln (k (n+1)) + 1"
by simp
also from Suc.IH have "(∑i = 1..n. ln (k i) + 1) = real (n+1) * ln (real (n+1))" ..
also from mvt[of "n+1"] have "… = real (n+2) * ln (real (n+2)) - ln (k (n+1)) - 1"
by simp
finally show ?case
by simp
qed simp
have **: "abs((∑i=1..n+1. ln i) - ((n+1) * ln (n+1) - (n+1))) ≤ 1 + ln(n+1)" for n::nat
proof -
have "(∑i=1..n+1. ln i) ≤ (∑i=1..n. ln i) + ln (n+1)"
by simp
also have "(∑i=1..n. ln i) ≤ (∑i=1..n. ln (k i))"
by (intro sum_mono, subst ln_le_cancel_iff) (auto simp: Suc_le_eq dest: lb ub)
also have "… = (∑i=1..n. ln (k i) + 1) - n"
also from * have "… = (n+1) * ln (n+1) - n"
by simp
finally have a_minus_b: "(∑i=1..n+1. ln i) - ((n+1) * ln (n+1) - (n+1)) ≤ 1 + ln (n+1)"
by simp

from * have "(n+1) * ln (n+1) - n = (∑i=1..n. ln (k i) + 1) - n"
by simp
also have "… = (∑i=1..n. ln (k i))"
also have "… ≤ (∑i=1..n. ln (i+1))"
by (intro sum_mono, subst ln_le_cancel_iff) (auto simp: Suc_le_eq dest: lb ub)
also from sum.shift_bounds_cl_nat_ivl[of "ln" 1 1 n] have "… = (∑i=1+1..n+1. ln i)" ..
also have "… = (∑i=1..n+1. ln i)"
by (rule sum.mono_neutral_left) auto
finally have b_minus_a: "((n+1) * ln (n+1) - (n+1)) - (∑i=1..n+1. ln i) ≤ 1"
by simp
have "0 ≤ ln (n+1)"
by simp
with b_minus_a have "((n+1) * ln (n+1) - (n+1)) - (∑i=1..n+1. ln i) ≤ 1 + ln (n+1)"
by linarith
with a_minus_b show ?thesis
by linarith
qed
from ‹n > 0› have "n ≥ 1" by simp
thus ?thesis
proof (induction n rule: dec_induct)
case base
then show ?case by simp
next
case (step n)
from ln_fact[of "n+1"] **[of n] show ?case by simp
qed
qed

lemma ln_fact_diff_bounds:
"abs(ln (fact n) - 2 * ln (fact (n div 2)) - n * ln 2) ≤ 4 * ln (if n = 0 then 1 else n) + 3"
proof (cases "n div 2 = 0")
case True
hence "n ≤ 1" by simp
with ln_le_minus_one[of "2::real"] show ?thesis by (cases n) simp_all
next
case False
then have "n > 1" by simp
let ?a = "real n * ln 2"
let ?b = "4 * ln (real n) + 3"
let ?l1 = "ln (fact (n div 2))"
let ?a1 = "real (n div 2) * ln (real (n div 2)) - real (n div 2)"
let ?b1 = "1 + ln (real (n div 2))"
let ?l2 = "ln (fact n)"
let ?a2 = "real n * ln (real n) - real n"
let ?b2 = "1 + ln (real n)"
have abs_a: "abs(?a - (?a2 - 2 * ?a1)) ≤ ?b - 2 * ?b1 - ?b2"
proof (cases "even n")
case True
then have "real (2 * (n div 2)) = real n"
by simp
then have n_div_2: "real (n div 2) = real n / 2"
by simp
from ‹n > 1› have *: "abs(?a - (?a2 - 2 * ?a1)) = 0"
by (simp add: n_div_2 ln_div algebra_simps)
from ‹even n› and ‹n > 1› have "0 ≤ ln (real n) - ln (real (n div 2))"
by (auto elim: evenE)
also have "2 * … ≤ 3 * ln (real n) - 2 * ln (real (n div 2))"
using ‹n > 1› by (auto intro!: ln_ge_zero)
also have "… = ?b - 2 * ?b1 - ?b2" by simp
finally show ?thesis using * by simp
next
case False
then have "real (2 * (n div 2)) = real (n - 1)"
by simp
with ‹n > 1› have n_div_2: "real (n div 2) = (real n - 1) / 2"
by simp
from ‹odd n› ‹n div 2 ≠ 0› have "n ≥ 3"
by presburger

have "?a - (?a2 - 2 * ?a1) = real n * ln 2 - real n * ln (real n) + real n +
2 * real (n div 2) * ln (real (n div 2)) - 2* real (n div 2)"
also from n_div_2 have "2 * real (n div 2) = real n - 1"
by simp
also have "real n * ln 2 - real n * ln (real n) + real n +
(real n - 1) * ln (real (n div 2)) - (real n - 1)
= real n * (ln (real n - 1) - ln (real n)) - ln (real (n div 2)) + 1"
using ‹n > 1› by (simp add: algebra_simps n_div_2 ln_div)
finally have lhs: "abs(?a - (?a2 - 2 * ?a1)) =
abs(real n * (ln (real n - 1) - ln (real n)) - ln (real (n div 2)) + 1)"
by simp

from ‹n > 1› have "real n * (ln (real n - 1) - ln (real n)) ≤ 0"
moreover from ‹n > 1› have "ln (real (n div 2)) ≤ ln (real n)" by simp
moreover {
have "exp 1 ≤ (3::real)" by (rule exp_le)
also from ‹n ≥ 3› have "… ≤ exp (ln (real n))" by simp
finally have "ln (real n) ≥ 1" by simp
}
ultimately have ub: "real n * (ln (real n - 1) - ln (real n)) - ln(real (n div 2)) + 1 ≤
3 * ln (real n) - 2 * ln(real (n div 2))" by simp

have mon: "real n' * (ln (real n') - ln (real n' - 1)) ≤
real n * (ln (real n) - ln (real n - 1))"
if "n ≥ 3" "n' ≥ n" for n n'::nat
proof (rule DERIV_nonpos_imp_nonincreasing[where f = "λx. x * (ln x - ln (x - 1))"])
fix t assume t: "real n ≤ t" "t ≤ real n'"
with that have "1 / (t - 1) ≥ ln (1 + 1/(t - 1))"
also from t that have "ln (1 + 1/(t - 1)) = ln t- ln (t - 1)"
by (simp add: ln_div [symmetric] field_simps)
finally have "ln t - ln (t - 1) ≤ 1 / (t - 1)" .
with that t
show "∃y. ((λx. x * (ln x - ln (x - 1))) has_field_derivative y) (at t) ∧ y ≤ 0"
by (intro exI[of _ "1 / (1 - t) + ln t - ln (t - 1)"])
(force intro!: derivative_eq_intros simp: field_simps)+
qed (use that in simp_all)

from ‹n > 1› have "ln 2 = ln (real n) - ln (real n / 2)"
also from ‹n > 1› have "… ≤ ln (real n) - ln (real (n div 2))"
by simp
finally have *: "3*ln 2 + ln(real (n div 2)) ≤ 3* ln(real n) - 2* ln(real (n div 2))"
by simp

have "- real n * (ln (real n - 1) - ln (real n)) + ln(real (n div 2)) - 1 =
real n * (ln (real n) - ln (real n - 1)) - 1 + ln(real (n div 2))"
also have "real n * (ln (real n) - ln (real n - 1)) ≤ 3 * (ln 3 - ln (3 - 1))"
using mon[OF _ ‹n ≥ 3›] by simp
also {
have "Some (Float 3 (-1)) = ub_ln 1 3" by code_simp
from ub_ln(1)[OF this] have "ln 3 ≤ (1.6 :: real)" by simp
also have "1.6 - 1 / 3 ≤ 2 * (2/3 :: real)" by simp
also have "2/3 ≤ ln (2 :: real)" by (rule ln_2_ge')
finally have "ln 3 - 1 / 3 ≤ 2 * ln (2 :: real)" by simp
}
hence "3 * (ln 3 - ln (3 - 1)) - 1 ≤ 3 * ln (2 :: real)" by simp
also note *
finally have "- real n * (ln (real n - 1) - ln (real n)) + ln(real (n div 2)) - 1 ≤
3 * ln (real n) - 2 * ln (real (n div 2))" by simp
hence lhs': "abs(real n * (ln (real n - 1) - ln (real n)) - ln(real (n div 2)) + 1) ≤
3 * ln (real n) - 2 * ln (real (n div 2))"
using ub by simp
have rhs: "?b - 2 * ?b1 - ?b2 = 3* ln (real n) - 2 * ln (real (n div 2))"
by simp
from ‹n > 1› have "ln (real (n div 2)) ≤ 3* ln (real n) - 2* ln (real (n div 2))"
by simp
with rhs lhs lhs' show ?thesis
by simp
qed
then have minus_a: "-?a ≤ ?b - 2 * ?b1 - ?b2 - (?a2 - 2 * ?a1)"
by simp
from abs_a have a: "?a ≤ ?b - 2 * ?b1 - ?b2 + ?a2 - 2 * ?a1"
by (simp)
from ln_fact_bounds[of "n div 2"] False have abs_l1: "abs(?l1 - ?a1) ≤ ?b1"
then have minus_l1: "?a1 - ?l1 ≤ ?b1"
by linarith
from abs_l1 have l1: "?l1 - ?a1 ≤ ?b1"
by linarith
from ln_fact_bounds[of n] False have abs_l2: "abs(?l2 - ?a2) ≤ ?b2"
then have l2: "?l2 - ?a2 ≤ ?b2"
by simp
from abs_l2 have minus_l2: "?a2 - ?l2 ≤ ?b2"
by simp
from minus_a minus_l1 l2 have "?l2 - 2 * ?l1 - ?a ≤ ?b"
by simp
moreover from a l1 minus_l2 have "- ?l2 + 2 * ?l1 + ?a ≤ ?b"
by simp
ultimately have "abs((?l2 - 2*?l1) - ?a) ≤ ?b"
by simp
then show ?thesis
by simp
qed

lemma ln_primefact:
assumes "n ≠ (0::nat)"
shows   "ln n = (∑d=1..n. if primepow d ∧ d dvd n then ln (aprimedivisor d) else 0)"
(is "?lhs = ?rhs")
proof -
have "?rhs = (∑d∈{x ∈ {1..n}. primepow x ∧ x dvd n}. ln (real (aprimedivisor d)))"
unfolding primepow_factors_def by (subst sum.inter_filter [symmetric]) simp_all
also have "{x ∈ {1..n}. primepow x ∧ x dvd n} = primepow_factors n"
using assms by (auto simp: primepow_factors_def dest: dvd_imp_le primepow_gt_Suc_0)
finally have *: "(∑d∈primepow_factors n. ln (real (aprimedivisor d))) = ?rhs" ..
from in_prime_factors_imp_prime prime_gt_0_nat
have pf_pos: "⋀p. p∈#prime_factorization n ⟹ p > 0"
by blast
from ln_msetprod[of "prime_factorization n", OF pf_pos] assms
have "ln n = (∑p∈#prime_factorization n. ln p)"
also from * sum_prime_factorization_conv_sum_primepow_factors[of n ln, OF assms(1)]
have "… = ?rhs" by simp
finally show ?thesis .
qed

context
begin

private lemma divisors:
fixes x d::nat
assumes "x ∈ {1..n}"
assumes "d dvd x"
shows "∃k∈{1..n div d}. x = d * k"
proof -
from assms have "x ≤ n"
by simp
then have ub: "x div d ≤ n div d"
by (simp add: div_le_mono ‹x ≤ n›)
from assms have "1 ≤ x div d" by (auto elim!: dvdE)
with ub have "x div d ∈ {1..n div d}"
by simp
with ‹d dvd x› show ?thesis by (auto intro!: bexI[of _ "x div d"])
qed

lemma ln_fact_conv_mangoldt: "ln (fact n) = (∑d=1..n. mangoldt d * floor (n / d))"
proof -
have *: "(∑da=1..n. if primepow da ∧ da dvd d then ln (aprimedivisor da) else 0) =
(∑(da::nat)=1..d. if primepow da ∧ da dvd d then ln (aprimedivisor da) else 0)"
if d: "d ∈ {1..n}" for d
by (rule sum.mono_neutral_right, insert d) (auto dest: dvd_imp_le)
have "(∑d=1..n. ∑da=1..d. if primepow da ∧
da dvd d then ln (aprimedivisor da) else 0) =
(∑d=1..n. ∑da=1..n. if primepow da ∧
da dvd d then ln (aprimedivisor da) else 0)"
by (rule sum.cong) (insert *, simp_all)
also have "… = (∑da=1..n. ∑d=1..n. if primepow da ∧
da dvd d then ln (aprimedivisor da) else 0)"
by (rule sum.swap)
also have "… = sum (λd. mangoldt d * floor (n/d)) {1..n}"
proof (rule sum.cong)
fix d assume d: "d ∈ {1..n}"
have "(∑da = 1..n. if primepow d ∧ d dvd da then ln (real (aprimedivisor d)) else 0) =
(∑da = 1..n. if d dvd da then mangoldt d else 0)"
by (intro sum.cong) (simp_all add: mangoldt_def)
also have "… = mangoldt d * real (card {x. x ∈ {1..n} ∧ d dvd x})"
by (subst sum.inter_filter [symmetric]) (simp_all add: algebra_simps)
also {
have "{x. x ∈ {1..n} ∧ d dvd x} = {x. ∃k ∈{1..n div d}. x=k*d}"
proof safe
fix x assume "x ∈ {1..n}" "d dvd x"
thus "∃k∈{1..n div d}. x = k * d" using divisors[of x n d] by auto
next
fix x k assume k: "k ∈ {1..n div d}"
from k have "k * d ≤ n div d * d" by (intro mult_right_mono) simp_all
also have "n div d * d ≤ n div d * d + n mod d" by (rule le_add1)
also have "… = n" by simp
finally have "k * d ≤ n" .
thus "k * d ∈ {1..n}" using d k by auto
qed auto
also have "… = (λk. k*d)  {1..n div d}"
by fast
also have "card … = card {1..n div d}"
by (rule card_image) (simp add: inj_on_def)
also have "… = n div d"
by simp
also have "... = ⌊n / d⌋"
finally have "real (card {x. x ∈ {1..n} ∧ d dvd x}) = real_of_int ⌊n / d⌋"
by force
}
finally show "(∑da = 1..n. if primepow d ∧ d dvd da then ln (real (aprimedivisor d)) else 0) =
mangoldt d * real_of_int ⌊real n / real d⌋" .
qed simp_all
finally have "(∑d=1..n. ∑da=1..d. if primepow da ∧
da dvd d then ln (aprimedivisor da) else 0) =
sum (λd. mangoldt d * floor (n/d)) {1..n}" .
with ln_primefact have "(∑d=1..n. ln d) =
(∑d=1..n. mangoldt d * floor (n/d))"
by simp
with ln_fact show ?thesis
by simp
qed

end

context
begin

private lemma div_2_mult_2_bds:
fixes n d :: nat
assumes "d > 0"
shows "0 ≤ ⌊n / d⌋ - 2 * ⌊(n div 2) / d⌋" "⌊n / d⌋ - 2 * ⌊(n div 2) / d⌋ ≤ 1"
proof -
have "⌊2::real⌋ * ⌊(n div 2) / d⌋ ≤ ⌊2 * ((n div 2) / d)⌋"
by (rule le_mult_floor) simp_all
also from assms have "… ≤ ⌊n / d⌋" by (intro floor_mono) (simp_all add: field_simps)
finally show "0 ≤ ⌊n / d⌋ - 2 * ⌊(n div 2) / d⌋" by (simp add: algebra_simps)
next
have "real (n div d) ≤ real (2 * ((n div 2) div d) + 1)"
by (subst div_mult2_eq [symmetric], simp only: mult.commute, subst div_mult2_eq) simp
thus "⌊n / d⌋ - 2 * ⌊(n div 2) / d⌋ ≤ 1"
unfolding of_nat_add of_nat_mult floor_conv_div_nat [symmetric] by simp_all
qed

private lemma n_div_d_eq_1: "d ∈ {n div 2 + 1..n} ⟹ ⌊real n / real d⌋ = 1"
by (cases "n = d") (auto simp: field_simps intro: floor_eq)

lemma psi_bounds_ln_fact:
shows "ln (fact n) - 2 * ln (fact (n div 2)) ≤ psi n"
"psi n - psi (n div 2) ≤ ln (fact n) - 2 * ln (fact (n div 2))"
proof -
fix n::nat
let ?k = "n div 2" and ?d = "n mod 2"
have *: "⌊?k / d⌋ = 0" if "d > ?k" for d
proof -
from that div_less have "0 = ?k div d" by simp
also have "… = ⌊?k / d⌋" by (rule floor_divide_of_nat_eq [symmetric])
finally show "⌊?k / d⌋ = 0" by simp
qed
have sum_eq: "(∑d=1..2*?k+?d. mangoldt d * ⌊?k / d⌋) = (∑d=1..?k. mangoldt d * ⌊?k / d⌋)"
by (intro sum.mono_neutral_right) (auto simp: *)
from ln_fact_conv_mangoldt have "ln (fact n) = (∑d=1..n. mangoldt d * ⌊n / d⌋)" .
also have "… = (∑d=1..n. mangoldt d * ⌊(2 * (n div 2) + n mod 2) / d⌋)"
by simp
also have "… ≤ (∑d=1..n. mangoldt d * (2 * ⌊?k / d⌋ + 1))"
using div_2_mult_2_bds(2)[of _ n]
by (intro sum_mono mult_left_mono, subst of_int_le_iff)
(auto simp: algebra_simps mangoldt_nonneg)
also have "… = 2 * (∑d=1..n. mangoldt d * ⌊(n div 2) / d⌋) + (∑d=1..n. mangoldt d)"
by (simp add: algebra_simps sum.distrib sum_distrib_left)
also have "… = 2 * (∑d=1..2*?k+?d. mangoldt d * ⌊(n div 2) / d⌋) + (∑d=1..n. mangoldt d)"
by presburger
also from sum_eq have "… = 2 * (∑d=1..?k. mangoldt d * ⌊(n div 2) / d⌋) + (∑d=1..n. mangoldt d)"
by presburger
also from ln_fact_conv_mangoldt psi_def have "… = 2 * ln (fact ?k) + psi n"
by presburger
finally show "ln (fact n) - 2 * ln (fact (n div 2)) ≤ psi n"
by simp
next
fix n::nat
let ?k = "n div 2" and  ?d = "n mod 2"
from psi_def have "psi n - psi ?k = (∑d=1..2*?k+?d. mangoldt d) - (∑d=1..?k. mangoldt d)"
by presburger
also have "… = sum mangoldt ({1..2 * (n div 2) + n mod 2} - {1..n div 2})"
by (subst sum_diff) simp_all
also have "… = (∑d∈({1..2 * (n div 2) + n mod 2} - {1..n div 2}).
(if d ≤ ?k then 0 else mangoldt d))"
by (intro sum.cong) simp_all
also have "… = (∑d=1..2*?k+?d. (if d ≤ ?k then 0 else mangoldt d))"
by (intro sum.mono_neutral_left) auto
also have "… = (∑d=1..n. (if d ≤ ?k then 0 else mangoldt d))"
by presburger
also have "… = (∑d=1..n. (if d ≤ ?k then mangoldt d * 0 else mangoldt d))"
by (intro sum.cong) simp_all
also from div_2_mult_2_bds(1) have "… ≤ (∑d=1..n. (if d ≤ ?k then mangoldt d * (⌊n/d⌋ - 2 * ⌊?k/d⌋) else mangoldt d))"
by (intro sum_mono)
(auto simp: algebra_simps mangoldt_nonneg intro!: mult_left_mono simp del: of_int_mult)
also from n_div_d_eq_1 have "… = (∑d=1..n. (if d ≤ ?k then mangoldt d * (⌊n/d⌋ - 2 * ⌊?k/d⌋) else mangoldt d * ⌊n/d⌋))"
by (intro sum.cong refl) auto
also have "… = (∑d=1..n. mangoldt d * real_of_int (⌊real n / real d⌋) -
(if d ≤ ?k then 2 * mangoldt d * real_of_int ⌊real ?k / real d⌋ else 0))"
by (intro sum.cong refl) (auto simp: algebra_simps)
also have "… = (∑d=1..n. mangoldt d * real_of_int (⌊real n / real d⌋)) -
(∑d=1..n. (if d ≤ ?k then 2 * mangoldt d * real_of_int ⌊real ?k / real d⌋ else 0))"
by (rule sum_subtractf)
also have "(∑d=1..n. (if d ≤ ?k then 2 * mangoldt d * real_of_int ⌊real ?k / real d⌋ else 0)) =
(∑d=1..?k. (if d ≤ ?k then 2 * mangoldt d * real_of_int ⌊real ?k / real d⌋ else 0))"
by (intro sum.mono_neutral_right) auto
also have "… = (∑d=1..?k. 2 * mangoldt d * real_of_int ⌊real ?k / real d⌋)"
by (intro sum.cong) simp_all
also have "… = 2 * (∑d=1..?k. mangoldt d * real_of_int ⌊real ?k / real d⌋)"
also have "(∑d = 1..n. mangoldt d * real_of_int ⌊real n / real d⌋) - … =
ln (fact n) - 2 * ln (fact (n div 2))"
finally show "psi n - psi (n div 2) ≤ ln (fact n) - 2 * ln (fact (n div 2))" .
qed

end

lemma psi_bounds_induct:
"real n * ln 2 - (4 * ln (real (if n = 0 then 1 else n)) + 3) ≤ psi n"
"psi n - psi (n div 2) ≤ real n * ln 2 + (4 * ln (real (if n = 0 then 1 else n)) + 3)"
proof -
from le_imp_neg_le[OF ln_fact_diff_bounds]
have "n * ln 2 - (4 * ln (if n = 0 then 1 else n) + 3)
≤ n * ln 2 - abs(ln (fact n) - 2 * ln (fact (n div 2)) - n * ln 2)"
by simp
also have "… ≤ ln (fact n) - 2 * ln (fact (n div 2))"
by simp
also from psi_bounds_ln_fact (1) have "… ≤ psi n"
by simp
finally show "real n * ln 2 - (4 * ln (real (if n = 0 then 1 else n)) + 3) ≤ psi n" .
next
from psi_bounds_ln_fact (2) have "psi n - psi (n div 2) ≤ ln (fact n) - 2 * ln (fact (n div 2))" .
also have "… ≤ n * ln 2 + abs(ln (fact n) - 2 * ln (fact (n div 2)) - n * ln 2)"
by simp
also from ln_fact_diff_bounds [of n]
have "abs(ln (fact n) - 2 * ln (fact (n div 2)) - n * ln 2)
≤ (4 * ln (real (if n = 0 then 1 else n)) + 3)" by simp
finally show "psi n - psi (n div 2) ≤ real n * ln 2 + (4 * ln (real (if n = 0 then 1 else n)) + 3)"
by simp
qed

subsection ‹Bounding the psi function›

text ‹
In this section, we will first prove the relatively tight estimate
@{prop "psi n ≤ 3 / 2 + ln 2 * n"} for @{term "n ≤ 128"} and then use the
recurrence we have just derived to extend it to @{prop "psi n ≤ 551 / 256"} for
@{term "n ≤ 1024"}, at which point applying the recurrence can be used to prove
the same bound for arbitrarily big numbers.

First of all, we will prove the bound for @{term "n <= 128"} using reflection and
approximation.
›

context
begin

private lemma Ball_insertD:
assumes "∀x∈insert y A. P x"
shows   "P y" "∀x∈A. P x"
using assms by auto

private lemma meta_eq_TrueE: "PROP A ≡ Trueprop True ⟹ PROP A"
by simp

private lemma pre_mangoldt_pos: "pre_mangoldt n > 0"
unfolding pre_mangoldt_def by (auto simp: primepow_gt_Suc_0)

private lemma psi_conv_pre_mangoldt: "psi n = ln (real (prod pre_mangoldt {1..n}))"
by (auto simp: psi_def mangoldt_def pre_mangoldt_def ln_prod primepow_gt_Suc_0 intro!: sum.cong)

private lemma eval_psi_aux1: "psi 0 = ln (real (numeral Num.One))"

private lemma eval_psi_aux2:
assumes "psi m = ln (real (numeral x))" "pre_mangoldt n = y" "m + 1 = n" "numeral x * y = z"
shows   "psi n = ln (real z)"
proof -
from assms(2) [symmetric] have [simp]: "y > 0" by (simp add: pre_mangoldt_pos)
have "psi n = psi (Suc m)" by (simp add: assms(3) [symmetric])
also have "… = ln (real y * (∏x = Suc 0..m. real (pre_mangoldt x)))"
using assms(2,3) [symmetric] by (simp add: psi_conv_pre_mangoldt prod.nat_ivl_Suc' mult_ac)
also have "… = ln (real y) + psi m"
by (subst ln_mult) (simp_all add: pre_mangoldt_pos prod_pos psi_conv_pre_mangoldt)
also have "psi m = ln (real (numeral x))" by fact
also have "ln (real y) + … = ln (real (numeral x * y))" by (simp add: ln_mult)
finally show ?thesis by (simp add: assms(4) [symmetric])
qed

private lemma Ball_atLeast0AtMost_doubleton:
assumes "psi 0 ≤ 3 / 2 * ln 2 * real 0"
assumes "psi 1 ≤ 3 / 2 * ln 2 * real 1"
shows   "(∀x∈{0..1}. psi x ≤ 3 / 2 * ln 2 * real x)"
using assms unfolding One_nat_def atLeast0_atMost_Suc ball_simps by auto

private lemma Ball_atLeast0AtMost_insert:
assumes "(∀x∈{0..m}. psi x ≤ 3 / 2 * ln 2 * real x)"
assumes "psi (numeral n) ≤ 3 / 2 * ln 2 * real (numeral n)" "m = pred_numeral n"
shows   "(∀x∈{0..numeral n}. psi x ≤ 3 / 2 * ln 2 * real x)"
using assms
by (subst numeral_eq_Suc[of n], subst atLeast0_atMost_Suc,
subst ball_simps, simp only: numeral_eq_Suc [symmetric])

private lemma eval_psi_ineq_aux:
assumes "psi n = x" "x ≤ 3 / 2 * ln 2 * n"
shows   "psi n ≤ 3 / 2 * ln 2 * n"
using assms by simp_all

private lemma eval_psi_ineq_aux2:
assumes "numeral m ^ 2 ≤ (2::nat) ^ (3 * n)"
shows   "ln (real (numeral m)) ≤ 3 / 2 * ln 2 * real n"
proof -
have "ln (real (numeral m)) ≤ 3 / 2 * ln 2 * real n ⟷
2 * log 2 (real (numeral m)) ≤ 3 * real n"
also have "2 * log 2 (real (numeral m)) = log 2 (real (numeral m ^ 2))"
by (subst of_nat_power, subst log_nat_power) simp_all
also have "… ≤ 3 * real n ⟷ real ((numeral m) ^ 2) ≤ 2 powr real (3 * n)"
by (subst Transcendental.log_le_iff) simp_all
also have "2 powr (3 * n) = real (2 ^ (3 * n))"
also have "real ((numeral m) ^ 2) ≤ … ⟷ numeral m ^ 2 ≤ (2::nat) ^ (3 * n)"
by (rule of_nat_le_iff)
finally show ?thesis using assms by blast
qed

private lemma eval_psi_ineq_aux_mono:
assumes "psi n = x" "psi m = x" "psi n ≤ 3 / 2 * ln 2 * n" "n ≤ m"
shows   "psi m ≤ 3 / 2 * ln 2 * m"
proof -
from assms have "psi m = psi n" by simp
also have "… ≤ 3 / 2 * ln 2 * n" by fact
also from ‹n ≤ m› have "… ≤ 3 / 2 * ln 2 * m" by simp
finally show ?thesis .
qed

lemma not_primepow_1_nat: "¬primepow (1 :: nat)" by auto

ML_file ‹bertrand.ML›

(* This should not take more than 1 minute *)
local_setup ‹fn ctxt =>
let
fun tac {context = ctxt, ...} =
let
val psi_cache = Bertrand.prove_psi ctxt 129
fun prove_psi_ineqs ctxt =
let
fun tac {context = ctxt, ...} =
HEADGOAL (resolve_tac ctxt @{thms eval_psi_ineq_aux2} THEN' Simplifier.simp_tac ctxt)
fun prove_by_approx n thm =
let
val thm = thm RS @{thm eval_psi_ineq_aux}
val [prem] = Thm.prems_of thm
val prem = Goal.prove ctxt [] [] prem tac
in
prem RS thm
end
fun prove_by_mono last_thm last_thm' thm =
let
val thm = @{thm eval_psi_ineq_aux_mono} OF [last_thm, thm, last_thm']
val [prem] = Thm.prems_of thm
val prem = Goal.prove ctxt [] [] prem (K (HEADGOAL (Simplifier.simp_tac ctxt)))
in
prem RS thm
end
fun go _ acc [] = acc
| go last acc ((n, x, thm) :: xs) =
let
val thm' =
case last of
NONE => prove_by_approx n thm
| SOME (last_x, last_thm, last_thm') =>
if last_x = x then
prove_by_mono last_thm last_thm' thm
else
prove_by_approx n thm
in
go (SOME (x, thm, thm')) (thm' :: acc) xs
end
in
rev o go NONE []
end

val psi_ineqs = prove_psi_ineqs ctxt psi_cache
fun prove_ball ctxt (thm1 :: thm2 :: thms) =
let
val thm = @{thm Ball_atLeast0AtMost_doubleton} OF [thm1, thm2]
fun solve_prem thm =
let
fun tac {context = ctxt, ...} = HEADGOAL (Simplifier.simp_tac ctxt)
val thm' = Goal.prove ctxt [] [] (Thm.cprem_of thm 1 |> Thm.term_of) tac
in
thm' RS thm
end
fun go thm thm' = (@{thm Ball_atLeast0AtMost_insert} OF [thm', thm]) |> solve_prem
in
fold go thms thm
end
| prove_ball _ _ = raise Match
in
HEADGOAL (resolve_tac ctxt [prove_ball ctxt psi_ineqs])
end
val thm = Goal.prove @{context} [] [] @{prop "∀n∈{0..128}. psi n ≤ 3 / 2 * ln 2 * n"} tac
in
Local_Theory.note ((@{binding psi_ubound_log_128}, []), [thm]) ctxt |> snd
end
›

end

context
begin

private lemma psi_ubound_aux:
defines "f ≡ λx::real. (4 * ln x + 3) / (ln 2 * x)"
assumes "x ≥ 2" "x ≤ y"
shows   "f x ≥ f y"
using assms(3)
proof (rule DERIV_nonpos_imp_nonincreasing, goal_cases)
case (1 t)
define f' where "f' = (λx. (1 - 4 * ln x) / x^2 / ln 2 :: real)"
from 1 assms(2) have "(f has_real_derivative f' t) (at t)" unfolding f_def f'_def
by (auto intro!: derivative_eq_intros simp: field_simps power2_eq_square)
moreover {
from ln_2_ge have "1/4 ≤ ln (2::real)" by simp
also from assms(2) 1 have "… ≤ ln t" by simp
finally have "ln t ≥ 1/4" .
}
with 1 assms(2) have "f' t ≤ 0" by (simp add: f'_def field_simps)
ultimately show ?case by (intro exI[of _ "f' t"]) simp_all
qed

text ‹
These next rules are used in combination with @{thm psi_bounds_induct} and
@{thm psi_ubound_log_128} to extend the upper bound for @{term "psi"} from values no greater
than 128 to values no greater than 1024. The constant factor of the upper bound changes every
time, but once we have reached 1024, the recurrence is self-sustaining in the sense that we do
not have to adjust the constant factor anymore in order to double the range.
›
lemma psi_ubound_log_double_cases':
assumes "⋀n. n ≤ m ⟹ psi n ≤ c * ln 2 * real n" "n ≤ m'" "m' = 2*m"
"c ≤ c'" "c ≥ 0" "m ≥ 1" "c' ≥ 1 + c/2 + (4 * ln (m+1) + 3) / (ln 2 * (m+1))"
shows   "psi n ≤ c' * ln 2 * real n"
proof (cases "n > m")
case False
hence "psi n ≤ c * ln 2 * real n" by (intro assms) simp_all
also have "c ≤ c'" by fact
finally show ?thesis by - (simp_all add: mult_right_mono)
next
case True
hence n: "n ≥ m+1" by simp
from psi_bounds_induct(2)[of n] True
have "psi n ≤ real n * ln 2 + 4 * ln (real n) + 3 + psi (n div 2)" by simp
also from assms have "psi (n div 2) ≤ c * ln 2 * real (n div 2)"
by (intro assms) simp_all
also have "real (n div 2) ≤ real n / 2" by simp
also have "c * ln 2 * … = c / 2 * ln 2 * real n" by simp
also have "real n * ln 2 + 4 * ln (real n) + 3 + … =
(1 + c/2) * ln 2 * real n + (4 * ln (real n) + 3)" by (simp add: field_simps)
also {
have "(4 * ln (real n) + 3) / (ln 2 * (real n)) ≤ (4 * ln (m+1) + 3) / (ln 2 * (m+1))"
using n assms by (intro psi_ubound_aux) simp_all
also from assms have "(4 * ln (m+1) + 3) / (ln 2 * (m+1)) ≤ c' - 1 - c/2"
finally have "4 * ln (real n) + 3 ≤ (c' - 1 - c/2) * ln 2 * real n"
using n by (simp add: field_simps)
}
also have "(1 + c / 2) * ln 2 * real n + (c' - 1 - c / 2) * ln 2 * real n = c' * ln 2 * real n"
finally show ?thesis using ‹c ≥ 0› by (simp_all add: mult_left_mono)
qed

end

lemma psi_ubound_log_double_cases:
assumes "∀n≤m. psi n ≤ c * ln 2 * real n"
"c' ≥ 1 + c/2 + (4 * ln (m+1) + 3) / (ln 2 * (m+1))"
"m' = 2*m" "c ≤ c'" "c ≥ 0" "m ≥ 1"
shows   "∀n≤m'. psi n ≤ c' * ln 2 * real n"
using assms(1) by (intro allI impI assms psi_ubound_log_double_cases'[of m c _ m' c']) auto

lemma psi_ubound_log_1024:
"∀n≤1024. psi n ≤ 551 / 256 * ln 2 * real n"
proof -
from psi_ubound_log_128 have "∀n≤128. psi n ≤ 3 / 2 * ln 2 * real n" by simp
hence "∀n≤256. psi n ≤ 1025 / 512 * ln 2 * real n"
proof (rule psi_ubound_log_double_cases, goal_cases)
case 1
have "Some (Float 624 (- 7)) = ub_ln 9 129" by code_simp
from ub_ln(1)[OF this] and ln_2_ge show ?case by (simp add: field_simps)
qed simp_all
hence "∀n≤512. psi n ≤ 549 / 256 * ln 2 * real n"
proof (rule psi_ubound_log_double_cases, goal_cases)
case 1
have "Some (Float 180 (- 5)) = ub_ln 7 257" by code_simp
from ub_ln(1)[OF this] and ln_2_ge show ?case by (simp add: field_simps)
qed simp_all
thus "∀n≤1024. psi n ≤ 551 / 256 * ln 2 * real n"
proof (rule psi_ubound_log_double_cases, goal_cases)
case 1
have "Some (Float 203 (- 5)) = ub_ln 7 513" by code_simp
from ub_ln(1)[OF this] and ln_2_ge show ?case by (simp add: field_simps)
qed simp_all
qed

lemma psi_bounds_sustained_induct:
assumes "4 * ln (1 + 2 ^ j) + 3 ≤ d * ln 2 * (1 + 2^j)"
assumes "4 / (1 + 2^j) ≤ d * ln 2"
assumes "0 ≤ c"
assumes "c / 2 + d + 1 ≤ c"
assumes "j ≤ k"
assumes "⋀n. n ≤ 2^k ⟹ psi n ≤ c * ln 2 * n"
assumes "n ≤ 2^(Suc k)"
shows "psi n ≤ c * ln 2 * n"
proof (cases "n ≤ 2^k")
case True
with assms(6) show ?thesis .
next
case False
from psi_bounds_induct(2)
have "psi n - psi (n div 2) ≤ real n * ln 2 + (4 * ln (real (if n = 0 then 1 else n)) + 3)" .
also from False have "(if n = 0 then 1 else n) = n"
by simp
finally have "psi n ≤ real n * ln 2 + (4 * ln (real n) + 3) + psi (n div 2)"
by simp
also from assms(6,7) have "psi (n div 2) ≤ c * ln 2 * (n div 2)"
by simp
also have "real (n div 2) ≤ real n / 2"
by simp
also have "real n * ln 2 + (4 * ln (real n) + 3) + c * ln 2 * (n / 2) ≤ c * ln 2 * real n"
proof (rule overpower_lemma[of
"λx. x * ln 2 + (4 * ln x + 3) + c * ln 2 * (x / 2)" "1+2^j"
"λx. c * ln 2 * x" "λx. c * ln 2 - ln 2 - 4 / x - c / 2 * ln 2"
"real n"])
from assms(1) have "4 * ln (1 + 2^j) + 3 ≤ d * ln 2 * (1 + 2^j)" .
also from assms(4) have "d ≤ c - c/2 - 1"
by simp
also have "(…) * ln 2 * (1 + 2 ^ j) = c * ln 2 * (1 + 2 ^ j) - c / 2 * ln 2 * (1 + 2 ^ j)
- (1 + 2 ^ j) * ln 2"
finally have "4 * ln (1 + 2^j) + 3 ≤ c * ln 2 * (1 + 2 ^ j) - c / 2 * ln 2 * (1 + 2 ^ j)
- (1 + 2 ^ j) * ln 2"
then show "(1 + 2 ^ j) * ln 2 + (4 * ln (1 + 2 ^ j) + 3)
+ c * ln 2 * ((1 + 2 ^ j) / 2) ≤ c * ln 2 * (1 + 2 ^ j)"
by simp
next
fix x::real
assume x: "1 + 2^j ≤ x"
moreover have "1 + 2 ^ j > (0::real)" by (simp add: add_pos_pos)
ultimately have x_pos: "x > 0" by linarith
show "((λx. c * ln 2 * x - (x * ln 2 + (4 * ln x + 3) + c * ln 2 * (x / 2)))
has_real_derivative c * ln 2 - ln 2 - 4 / x - c / 2 * ln 2) (at x)"
by (rule derivative_eq_intros refl | simp add: ‹0 < x›)+
from ‹0 < x› ‹0 < 1 + 2^j› have "0 < x * (1 + 2^j)"
by (rule mult_pos_pos)
have "4 / x ≤ 4 / (1 + 2^j)"
by (intro divide_left_mono mult_pos_pos add_pos_pos x x_pos) simp_all
also from assms(2) have "4 / (1 + 2^j) ≤ d * ln 2" .
also from assms(4) have "d ≤ c - c/2 - 1" by simp
also have "… * ln 2 = c * ln 2 - c/2 * ln 2 - ln 2" by (simp add: algebra_simps)
finally show "0 ≤ c * ln 2 - ln 2 - 4 / x - c / 2 * ln 2" by simp
next
have "1 + 2^j = real (1 + 2^j)" by simp
also from assms(5) have "… ≤ real (1 + 2^k)" by simp
also from False have "2^k ≤ n - 1" by simp
finally show "1 + 2^j ≤ real n" using False by simp
qed
finally show ?thesis using assms by - (simp_all add: mult_left_mono)
qed

lemma psi_bounds_sustained:
assumes "⋀n. n ≤ 2^k ⟹ psi n ≤ c * ln 2 * n"
assumes "4 * ln (1 + 2^k) + 3 ≤ (c/2 - 1) * ln 2 * (1 + 2^k)"
assumes "4 / (1 + 2^k) ≤ (c/2 - 1) * ln 2"
assumes "c ≥ 0"
shows "psi n ≤ c * ln 2 * n"
proof -
have *: "psi n ≤ c * ln 2 * n" if "n ≤ 2^j" for j n
using that
proof (induction j arbitrary: n)
case 0
with assms(4) 0 show ?case unfolding psi_def mangoldt_def by (cases n) auto
next
case (Suc j)
show ?case
proof (cases "k ≤ j")
case True
from assms(4) have c_div_2: "c/2 + (c/2 - 1) + 1 ≤ c"
by simp
from psi_bounds_sustained_induct[of k "c/2 -1" c j,
OF assms(2) assms(3) assms(4) c_div_2 True Suc.IH Suc.prems]
show ?thesis by simp
next
case False
then have j_lt_k: "Suc j ≤ k" by simp
from Suc.prems have "n ≤ 2 ^ Suc j" .
also have "(2::nat) ^ Suc j ≤ 2 ^ k"
using power_increasing[of "Suc j" k "2::nat", OF j_lt_k]
by simp
finally show ?thesis using assms(1) by simp
qed
qed
have "n < 2 ^ n" by (induction n) simp_all
with *[of n n] show ?thesis by simp
qed

lemma psi_ubound_log: "psi n ≤ 551 / 256 * ln 2 * n"
proof (rule psi_bounds_sustained)
show "0 ≤ 551 / (256 :: real)" by simp
next
fix n :: nat assume "n ≤ 2 ^ 10"
with psi_ubound_log_1024 show "psi n ≤ 551 / 256 * ln 2 * real n" by auto
next
have "4 / (1 + 2 ^ 10) ≤ (551 / 256 / 2 - 1) * (2/3 :: real)"
by simp
also have "… ≤ (551 / 256 / 2 - 1) * ln 2"
by (intro mult_left_mono ln_2_ge') simp_all
finally show "4 / (1 + 2 ^ 10) ≤ (551 / 256 / 2 - 1) * ln (2 :: real)" .
next
have "Some (Float 16 (-1)) = ub_ln 3 1025" by code_simp
from ub_ln(1)[OF this] and ln_2_ge
have "2048 * ln 1025 + 1536 ≤ 39975 * (ln 2::real)" by simp
thus "4 * ln (1 + 2 ^ 10) + 3 ≤ (551 / 256 / 2 - 1) * ln 2 * (1 + 2 ^ 10 :: real)"
by simp
qed

lemma psi_ubound_3_2: "psi n ≤ 3/2 * n"
proof -
have "(551 / 256) * ln 2 ≤ (551 / 256) * (16/23 :: real)"
by (intro mult_left_mono ln_2_le') auto
also have "… ≤ 3 / 2" by simp
finally have "551 / 256 * ln 2 ≤ 3/(2::real)" .
with of_nat_0_le_iff mult_right_mono have "551 / 256 * ln 2 * n ≤ 3/2 * n"
by blast
with psi_ubound_log[of "n"] show ?thesis
by linarith
qed

subsection ‹Doubling psi and theta›

lemma psi_residues_compare_2:
"psi_odd_2 n ≤ psi_even_2 n"
proof -
have "psi_odd_2 n = (∑d∈{d. d ∈ {2..n} ∧ primepow_odd d}. mangoldt_odd d)"
unfolding mangoldt_odd_def by (rule sum.mono_neutral_right) auto
also have "… = (∑d∈{d. d ∈ {2..n} ∧ primepow_odd d}. ln (real (aprimedivisor d)))"
by (intro sum.cong refl) (simp add: mangoldt_odd_def)
also have "… ≤ (∑d∈{d. d ∈ {2..n} ∧ primepow_even d}. ln (real (aprimedivisor d)))"
proof (rule sum_le_included [where i = "λy. y * aprimedivisor y"]; clarify?)
fix d :: nat assume "d ∈ {2..n}" "primepow_odd d"
note d = this
then obtain p k where d': "k ≥ 1" "prime p" "d = p ^ (2*k+1)"
by (auto simp: primepow_odd_def)
from d' have "p ^ (2 * k) ≤ p ^ (2 * k + 1)"
by (subst power_increasing_iff) (auto simp: prime_gt_Suc_0_nat)
also from d d' have "… ≤ n" by simp
finally have "p ^ (2 * k) ≤ n" .
moreover from d' have "p ^ (2 * k) > 1"
by (intro one_less_power) (simp_all add: prime_gt_Suc_0_nat)
ultimately have "p ^ (2 * k) ∈ {2..n}" by simp
moreover from d' have "primepow_even (p ^ (2 * k))"
by (auto simp: primepow_even_def)
ultimately show "∃y∈{d ∈ {2..n}. primepow_even d}. y * aprimedivisor y = d ∧
ln (real (aprimedivisor d)) ≤ ln (real (aprimedivisor y))" using d'
by (intro bexI[of _ "p ^ (2 * k)"])
(auto simp: aprimedivisor_prime_power aprimedivisor_primepow)
also have "… = (∑d∈{d. d ∈ {2..n} ∧ primepow_even d}. mangoldt_even d)"
by (intro sum.cong refl) (simp add: mangoldt_even_def)
also have "… = psi_even_2 n"
unfolding mangoldt_even_def by (rule sum.mono_neutral_left) auto
finally show ?thesis .
qed

lemma psi_residues_compare:
"psi_odd n ≤ psi_even n"
proof -
have "¬ primepow_odd 1" by (simp add: primepow_odd_def)
hence *: "mangoldt_odd 1 = 0" by (simp add: mangoldt_odd_def)
have "¬ primepow_even 1"
using primepow_gt_Suc_0[OF primepow_even_imp_primepow, of 1] by auto
with mangoldt_even_def have **: "mangoldt_even 1 = 0"
by simp
from psi_odd_def have "psi_odd n = (∑d=1..n. mangoldt_odd d)"
by simp
also from * have "… = psi_odd_2 n"
by (cases "n ≥ 1") (simp_all add: eval_nat_numeral sum.atLeast_Suc_atMost)
also from psi_residues_compare_2 have "… ≤ psi_even_2 n" .
also from ** have "… = psi_even n"
by (cases "n ≥ 1") (simp_all add: eval_nat_numeral sum.atLeast_Suc_atMost psi_even_def)
finally show ?thesis .
qed

lemma primepow_iff_even_sqr:
"primepow n ⟷ primepow_even (n^2)"
by (cases "n = 0")
(auto simp: primepow_even_altdef aprimedivisor_primepow_power primepow_power_iff_nat
prime_elem_multiplicity_power_distrib prime_aprimedivisor' prime_imp_prime_elem
unit_factor_nat_def primepow_gt_0_nat dest: primepow_gt_Suc_0)

lemma psi_sqrt: "psi (Discrete.sqrt n) = psi_even n"
proof (induction n)
case 0
with psi_def psi_even_def show ?case by simp
next
case (Suc n)
then show ?case
proof cases
assume asm: "∃ m. Suc n = m^2"
with sqrt_Suc have sqrt_seq: "Discrete.sqrt(Suc n) = Suc(Discrete.sqrt n)"
by simp
from asm obtain "m" where "   Suc n = m^2"
by blast
with sqrt_seq have "Suc(Discrete.sqrt n) = m"
by simp
with ‹Suc n = m^2› have suc_sqrt_n_sqrt: "(Suc(Discrete.sqrt n))^2 = Suc n"
by simp
from sqrt_seq have "psi (Discrete.sqrt (Suc n)) = psi (Suc (Discrete.sqrt n))"
by simp
also from psi_def have "… = psi (Discrete.sqrt n) + mangoldt (Suc (Discrete.sqrt n))"
by simp
also from Suc.IH have "psi (Discrete.sqrt n) = psi_even n" .
also have "mangoldt (Suc (Discrete.sqrt n)) = mangoldt_even (Suc n)"
proof (cases "primepow (Suc(Discrete.sqrt n))")
case True
with primepow_iff_even_sqr have True2: "primepow_even ((Suc(Discrete.sqrt n))^2)"
by simp
from suc_sqrt_n_sqrt have "mangoldt_even (Suc n) = mangoldt_even ((Suc(Discrete.sqrt n))^2)"
by simp
also from mangoldt_even_def True2
have "… = ln (aprimedivisor ((Suc (Discrete.sqrt n))^2))"
by simp
also from True have "aprimedivisor ((Suc (Discrete.sqrt n))^2) = aprimedivisor (Suc (Discrete.sqrt n))"
also from True have "ln (…) = mangoldt (Suc (Discrete.sqrt n))"
finally show ?thesis ..
next
case False
with primepow_iff_even_sqr
have False2: "¬ primepow_even ((Suc(Discrete.sqrt n))^2)"
by simp
from suc_sqrt_n_sqrt have "mangoldt_even (Suc n) = mangoldt_even ((Suc(Discrete.sqrt n))^2)"
by simp
also from mangoldt_even_def False2
have "… = 0"
by simp
also from False have "… = mangoldt (Suc (Discrete.sqrt n))"
finally show ?thesis ..
qed
also from psi_even_def have "psi_even n + mangoldt_even (Suc n) = psi_even (Suc n)"
by simp
finally show ?case .
next
assume asm: "¬(∃m. Suc n = m^2)"
with sqrt_Suc have sqrt_eq: "Discrete.sqrt (Suc n) = Discrete.sqrt n"
by simp
then have lhs: "psi (Discrete.sqrt (Suc n)) = psi (Discrete.sqrt n)"
by simp
have "¬ primepow_even (Suc n)"
proof
assume "primepow_even (Suc n)"
with primepow_even_def obtain "p" "k"
where "1 ≤ k ∧ prime p ∧ Suc n = p ^ (2 * k)"
by blast
with power_even_eq have "Suc n = (p ^ k)^2"
by simp
with asm show False by blast
qed
with psi_even_def mangoldt_even_def
have rhs: "psi_even (Suc n) = psi_even n"
by simp
from Suc.IH lhs rhs show ?case
by simp
qed
qed

lemma mangoldt_split:
"mangoldt d = mangoldt_1 d + mangoldt_even d + mangoldt_odd d"
proof (cases "primepow d")
case False
thus ?thesis
by (auto simp: mangoldt_def mangoldt_1_def mangoldt_even_def mangoldt_odd_def
dest: primepow_even_imp_primepow primepow_odd_imp_primepow)
next
case True
thus ?thesis
by (auto simp: mangoldt_def mangoldt_1_def mangoldt_even_def mangoldt_odd_def primepow_cases)
qed

lemma psi_split: "psi n = theta n + psi_even n + psi_odd n"
by (induction n)
(simp_all add: psi_def theta_def psi_even_def psi_odd_def mangoldt_1_def mangoldt_split)

lemma psi_mono: "m ≤ n ⟹ psi m ≤ psi n" unfolding psi_def
by (intro sum_mono2 mangoldt_nonneg) auto

lemma psi_pos: "0 ≤ psi n"
by (auto simp: psi_def intro!: sum_nonneg mangoldt_nonneg)

lemma mangoldt_odd_pos: "0 ≤ mangoldt_odd d"
using aprimedivisor_gt_Suc_0[of d]
by (auto simp: mangoldt_odd_def of_nat_le_iff[of 1, unfolded of_nat_1] Suc_le_eq
intro!: ln_ge_zero dest!: primepow_odd_imp_primepow primepow_gt_Suc_0)

lemma psi_odd_mono: "m ≤ n ⟹ psi_odd m ≤ psi_odd n"
using mangoldt_odd_pos sum_mono2[of "{1..n}" "{1..m}" "mangoldt_odd"]

lemma psi_odd_pos: "0 ≤ psi_odd n"
by (auto simp: psi_odd_def intro!: sum_nonneg mangoldt_odd_pos)

lemma psi_theta:
"theta n + psi (Discrete.sqrt n) ≤ psi n" "psi n ≤ theta n + 2 * psi (Discrete.sqrt n)"
using psi_odd_pos[of n] psi_residues_compare[of n] psi_sqrt[of n] psi_split[of n]
by simp_all

context
begin

private lemma sum_minus_one:
"(∑x ∈ {1..y}. (- 1 :: real) ^ (x + 1)) = (if odd y then 1 else 0)"
by (induction y) simp_all

private lemma div_invert:
fixes x y n :: nat
assumes "x > 0" "y > 0" "y ≤ n div x"
shows "x ≤ n div y"
proof -
from assms(1,3) have "y * x ≤ (n div x) * x"
by simp
also have "… ≤ n"
finally have "y * x ≤ n" .
with assms(2) show ?thesis
using div_le_mono[of "y*x" n y] by simp
qed

lemma sum_expand_lemma:
"(∑d=1..n. (-1) ^ (d + 1) * psi (n div d)) =
(∑d = 1..n. (if odd (n div d) then 1 else 0) * mangoldt d)"
proof -
have **: "x ≤ n" if "x ≤ n div y" for x y
using div_le_dividend order_trans that by blast
have "(∑d=1..n. (-1)^(d+1) * psi (n div d)) =
(∑d=1..n. (-1)^(d+1) * (∑e=1..n div d. mangoldt e))"
also have "… = (∑d = 1..n. ∑e = 1..n div d. (-1)^(d+1) * mangoldt e)"
also from ** have "… = (∑d = 1..n. ∑e∈{y∈{1..n}. y ≤ n div d}. (-1)^(d+1) * mangoldt e)"
by (intro sum.cong) auto
also have "… = (∑y = 1..n. ∑x | x ∈ {1..n} ∧ y ≤ n div x. (- 1) ^ (x + 1) * mangoldt y)"
by (rule sum.swap_restrict) simp_all
also have "… = (∑y = 1..n. ∑x | x ∈ {1..n} ∧ x ≤ n div y. (- 1) ^ (x + 1) * mangoldt y)"
by (intro sum.cong) (auto intro: div_invert)
also from ** have "… = (∑y = 1..n. ∑x ∈ {1..n div y}. (- 1) ^ (x + 1) * mangoldt y)"
by (intro sum.cong) auto
also have "… = (∑y = 1..n. (∑x ∈ {1..n div y}. (- 1) ^ (x + 1)) * mangoldt y)"
by (intro sum.cong) (simp_all add: sum_distrib_right)
also have "… = (∑y = 1..n. (if odd (n div y) then 1 else 0) * mangoldt y)"
by (intro sum.cong refl) (simp_all only: sum_minus_one)
finally show ?thesis .
qed

private lemma floor_half_interval:
fixes n d :: nat
assumes "d ≠ 0"
shows "real (n div d) - real (2 * ((n div 2) div d)) = (if odd (n div d) then 1 else 0)"
proof -
have "((n div 2) div d) = (n div (2 * d))"
by (rule div_mult2_eq[symmetric])
also have "… = ((n div d) div 2)"
also have "real (n div d) - real (2 * …) = (if odd (n div d) then 1 else 0)"
by (cases "odd (n div d)", cases "n div d = 0 ", simp_all)
finally show ?thesis by simp
qed

lemma fact_expand_psi:
"ln (fact n) - 2 * ln (fact (n div 2)) = (∑d=1..n. (-1)^(d+1) * psi (n div d))"
proof -
have "ln (fact n) - 2 * ln (fact (n div 2)) =
(∑d=1..n. mangoldt d * ⌊n / d⌋) - 2 * (∑d=1..n div 2. mangoldt d * ⌊(n div 2) / d⌋)"
also have "(∑d=1..n div 2. mangoldt d * ⌊real (n div 2) / d⌋) =
(∑d=1..n. mangoldt d * ⌊real (n div 2) / d⌋)"
by (rule sum.mono_neutral_left) (auto simp: floor_unique[of 0])
also have "2 * … = (∑d=1..n. mangoldt d * 2 * ⌊real (n div 2) / d⌋)"
also have "(∑d=1..n. mangoldt d * ⌊n / d⌋) - … =
(∑d=1..n. (mangoldt d * ⌊n / d⌋ - mangoldt d * 2 * ⌊real (n div 2) / d⌋))"
also have "… = (∑d=1..n. mangoldt d * (⌊n / d⌋ - 2 * ⌊real (n div 2) / d⌋))"
also have "… = (∑d=1..n. mangoldt d * (if odd(n div d) then 1 else 0))"
by (intro sum.cong refl)
(simp_all add: floor_conv_div_nat [symmetric] floor_half_interval [symmetric])
also have "… = (∑d=1..n. (if odd(n div d) then 1 else 0) * mangoldt d)"
also from sum_expand_lemma[symmetric] have "… = (∑d=1..n. (-1)^(d+1) * psi (n div d))" .
finally show ?thesis .
qed

end

lemma psi_expansion_cutoff:
assumes "m ≤ p"
shows   "(∑d=1..2*m. (-1)^(d+1) * psi (n div d)) ≤ (∑d=1..2*p. (-1)^(d+1) * psi (n div d))"
"(∑d=1..2*p+1. (-1)^(d+1) * psi (n div d)) ≤ (∑d=1..2*m+1. (-1)^(d+1) * psi (n div d))"
using assms
proof (induction m rule: inc_induct)
case (step k)
have "(∑d = 1..2 * k. (-1)^(d + 1) * psi (n div d)) ≤
(∑d = 1..2 * Suc k. (-1)^(d + 1) * psi (n div d))"
with step.IH(1)
show "(∑d = 1..2 * k. (-1)^(d + 1) * psi (n div d))
≤ (∑d = 1..2 * p. (-1)^(d + 1) * psi (n div d))"
by simp
from step.IH(2)
have "(∑d = 1..2 * p + 1. (-1)^(d + 1) * psi (n div d))
≤ (∑d = 1..2 * Suc k + 1. (-1)^(d + 1) * psi (n div d))" .
also have "… ≤ (∑d = 1..2 * k + 1. (-1)^(d + 1) * psi (n div d))"
finally show "(∑d = 1..2 * p + 1. (-1)^(d + 1) * psi (n div d))
≤ (∑d = 1..2 * k + 1. (-1)^(d + 1) * psi (n div d))" .
qed simp_all

lemma fact_psi_bound_even:
assumes "even k"
shows   "(∑d=1..k. (-1)^(d+1) * psi (n div d)) ≤ ln (fact n) - 2 * ln (fact (n div 2))"
proof -
have "(∑d=1..k. (-1)^(d+1) * psi (n div d)) ≤ (∑d = 1..n. (- 1) ^ (d + 1) * psi (n div d))"
proof (cases "k ≤ n")
case True
with psi_expansion_cutoff(1)[of "k div 2" "n div 2" n]
have "(∑d=1..2*(k div 2). (-1)^(d+1) * psi (n div d))
≤ (∑d = 1..2*(n div 2). (- 1) ^ (d + 1) * psi (n div d))"
by simp
also from assms have "2*(k div 2) = k"
by simp
also have "(∑d = 1..2*(n div 2). (- 1) ^ (d + 1) * psi (n div d))
≤ (∑d = 1..n. (- 1) ^ (d + 1) * psi (n div d))"
proof (cases "even n")
case True
then show ?thesis
by simp
next
case False
from psi_pos have "(∑d = 1..2*(n div 2). (- 1) ^ (d + 1) * psi (n div d))
≤ (∑d = 1..2*(n div 2) + 1. (- 1) ^ (d + 1) * psi (n div d))"
by simp
with False show ?thesis
by simp
qed
finally show ?thesis .
next
case False
hence *: "n div 2 ≤ (k-1) div 2"
by simp
have "(∑d=1..k. (-1)^(d+1) * psi (n div d)) ≤
(∑d=1..2*((k-1) div 2) + 1. (-1)^(d+1) * psi (n div d))"
proof (cases "k = 0")
case True
with psi_pos show ?thesis by simp
next
case False
with sum.cl_ivl_Suc[of "λd. (-1)^(d+1) * psi (n div d)" 1 "k-1"]
have "(∑d=1..k. (-1)^(d+1) * psi (n div d)) = (∑d=1..k-1. (-1)^(d+1) * psi (n div d))
+ (-1)^(k+1) * psi (n div k)"
by simp
also from assms psi_pos have "(-1)^(k+1) * psi (n div k) ≤ 0"
by simp
also from assms False have "k-1 = 2*((k-1) div 2) + 1"
by presburger
finally show ?thesis by simp
qed
also from * psi_expansion_cutoff(2)[of "n div 2" "(k-1) div 2" n]
have "… ≤ (∑d=1..2*(n div 2) + 1. (-1)^(d+1) * psi (n div d))" by blast
also have "… ≤ (∑d = 1..n. (- 1) ^ (d + 1) * psi (n div d))"
by (cases "even n") (simp_all add: psi_def)
finally show ?thesis .
qed
also from fact_expand_psi have "… = ln (fact n) - 2 * ln (fact (n div 2))" ..
finally show ?thesis .
qed

lemma fact_psi_bound_odd:
assumes "odd k"
shows "ln (fact n) - 2 * ln (fact (n div 2)) ≤ (∑d=1..k. (-1)^(d+1) * psi (n div d))"
proof -
from fact_expand_psi
have "ln (fact n) - 2 * ln (fact (n div 2)) = (∑d = 1..n. (- 1) ^ (d + 1) * psi (n div d))" .
also have "… ≤ (∑d=1..k. (-1)^(d+1) * psi (n div d))"
proof (cases "k ≤ n")
case True
have "(∑d=1..n. (-1)^(d+1) * psi (n div d)) ≤ (
∑d=1..2*(n div 2)+1. (-1)^(d+1) * psi (n div d))"
by (cases "even n") (simp_all add: psi_pos)
also from True assms psi_expansion_cutoff(2)[of "k div 2" "n div 2" n]
have "… ≤ (∑d=1..k. (-1)^(d+1) * psi (n div d))"
by simp
finally show ?thesis .
next
case False
have "(∑d=1..n. (-1)^(d+1) * psi (n div d)) ≤ (∑d=1..2*((n+1) div 2). (-1)^(d+1) * psi (n div d))"
by (cases "even n") (simp_all add: psi_def)
also from False assms psi_expansion_cutoff(1)[of "(n+1) div 2" "k div 2" n]
have "(∑d=1..2*((n+1) div 2). (-1)^(d+1) * psi (n div d)) ≤ (∑d=1..2*(k div 2). (-1)^(d+1) * psi (n div d))"
by simp
also from assms have "… ≤ (∑d=1..k. (-1)^(d+1) * psi (n div d))"
by (auto elim: oddE simp: psi_pos)
finally show ?thesis .
qed
finally show ?thesis .
qed

lemma fact_psi_bound_2_3:
"psi n - psi (n div 2) ≤ ln (fact n) - 2 * ln (fact (n div 2))"
"ln (fact n) - 2 * ln (fact (n div 2)) ≤ psi n - psi (n div 2) + psi (n div 3)"
proof -
show "psi n - psi (n div 2) ≤ ln (fact n) - 2 * ln (fact (n div 2))"
by (rule psi_bounds_ln_fact (2))
next
from fact_psi_bound_odd[of 3 n] have "ln (fact n) - 2 * ln (fact (n div 2))
≤ (∑d = 1..3. (- 1) ^ (d + 1) * psi (n div d))"
by simp
also have "… = psi n - psi (n div 2) + psi (n div 3)"
finally show "ln (fact n) - 2 * ln (fact (n div 2)) ≤ psi n - psi (n div 2) + psi (n div 3)" .
qed

lemma ub_ln_1200: "ln 1200 ≤ 57 / (8 :: real)"
proof -
have "Some (Float 57 (-3)) = ub_ln 8 1200" by code_simp
from ub_ln(1)[OF this] show ?thesis by simp
qed

lemma psi_double_lemma:
assumes "n ≥ 1200"
shows "real n / 6 ≤ psi n - psi (n div 2)"
proof -
from ln_fact_diff_bounds
have "¦ln (fact n) - 2 * ln (fact (n div 2)) - real n * ln 2¦
≤ 4 * ln (real (if n = 0 then 1 else n)) + 3" .
with assms have "ln (fact n) - 2 * ln (fact (n div 2))
≥ real n * ln 2 - 4 * ln (real n) - 3"
by simp
moreover have "real n * ln 2 - 4 * ln (real n) - 3 ≥ 2 / 3 * n"
proof (rule overpower_lemma[of "λn. 2/3 * n" 1200])
show "2 / 3 * 1200 ≤ 1200 * ln 2 - 4 * ln 1200 - (3::real)"
using ub_ln_1200 ln_2_ge by linarith
next
fix x::real
assume "1200 ≤ x"
then have "0 < x"
by simp
show "((λx. x * ln 2 - 4 * ln x - 3 - 2 / 3 * x)
has_real_derivative ln 2 - 4 / x - 2 / 3) (at x)"
by (rule derivative_eq_intros refl | simp add: ‹0 < x›)+
next
fix x::real
assume "1200 ≤ x"
then have "12 / x ≤ 12 / 1200" by simp
then have "0 ≤ 0.67 - 4 / x - 2 / 3" by simp
also have "0.67 ≤ ln (2::real)" using ln_2_ge by simp
finally show "0 ≤ ln 2 - 4 / x - 2 / 3" by simp
next
from assms show "1200 ≤ real n"
by simp
qed
ultimately have "2 / 3 * real n ≤ ln (fact n) - 2 * ln (fact (n div 2))"
by simp
with psi_ubound_3_2[of "n div 3"]
have "n/6 + psi (n div 3) ≤ ln (fact n) - 2 * ln (fact (n div 2))"
by simp
with fact_psi_bound_2_3[of "n"] show ?thesis
by simp
qed

lemma theta_double_lemma:
assumes "n ≥ 1200"
shows "theta (n div 2) < theta n"
proof -
from psi_theta[of "n div 2"] psi_pos[of "Discrete.sqrt (n div 2)"]
have theta_le_psi_n_2: "theta (n div 2) ≤ psi (n div 2)"
by simp
have "(Discrete.sqrt n * 18)^2 ≤ 324 * n"
by simp
from mult_less_cancel2[of "324" "n" "n"] assms have "324 * n < n^2"
with ‹(Discrete.sqrt n * 18)^2 ≤ 324 * n› have "(Discrete.sqrt n*18)^2 < n^2"
by presburger
with power2_less_imp_less assms have "Discrete.sqrt n * 18 < n"
by blast
with psi_ubound_3_2[of "Discrete.sqrt n"] have "2 * psi (Discrete.sqrt n) < n / 6"
by simp
with psi_theta[of "n"] have psi_lt_theta_n: "psi n - n / 6 < theta n"
by simp
from psi_double_lemma[OF assms(1)] have "psi (n div 2) ≤ psi n - n / 6"
by simp
with theta_le_psi_n_2 psi_lt_theta_n show ?thesis
by simp
qed

subsection ‹Proof of the main result›

lemma theta_mono: "mono theta"
by (auto simp: theta_def [abs_def] intro!: monoI sum_mono2)

lemma theta_lessE:
assumes "theta m < theta n" "m ≥ 1"
obtains p where "p ∈ {m<..n}" "prime p"
proof -
from mono_invE[OF theta_mono assms(1)] have "m ≤ n" by blast
hence "theta n = theta m + (∑p∈{m<..n}. if prime p then ln (real p) else 0)"
unfolding theta_def using assms(2)
by (subst sum.union_disjoint [symmetric]) (auto simp: ivl_disj_un)
also note assms(1)
finally have "(∑p∈{m<..n}. if prime p then ln (real p) else 0) ≠ 0" by simp
from sum.not_neutral_contains_not_neutral [OF this] guess p .
thus ?thesis using that[of p] by (auto intro!: exI[of _ p] split: if_splits)
qed

theorem bertrand:
fixes   n :: nat
assumes "n > 1"
shows   "∃p∈{n<..<2*n}. prime p"
proof cases
assume n_less: "n < 600"
define prime_constants
where "prime_constants = {2, 3, 5, 7, 13, 23, 43, 83, 163, 317, 631::nat}"
from ‹n > 1› n_less have "∃p ∈ prime_constants. n < p ∧ p < 2 * n"
unfolding bex_simps greaterThanLessThan_iff prime_constants_def by presburger
moreover have "∀p∈prime_constants. prime p"
unfolding prime_constants_def ball_simps HOL.simp_thms by (intro conjI; pratt (silent))
ultimately show ?thesis
unfolding greaterThanLessThan_def greaterThan_def lessThan_def by blast
next
assume n: "¬(n < 600)"
from n have "theta n < theta (2 * n)" using theta_double_lemma[of "2 * n"] by simp
with assms obtain p where "p ∈ {n<..2*n}" "prime p" by (auto elim!: theta_lessE)
moreover from assms have "¬prime (2*n)" by (auto dest!: prime_product)
with ‹prime p› have "p ≠ 2 * n" by auto
ultimately show ?thesis
by auto
qed

subsection ‹Proof of Mertens' first theorem›

text ‹
The following proof of Mertens' first theorem was ported from John Harrison's HOL Light
proof by Larry Paulson:
›

lemma sum_integral_ubound_decreasing':
fixes f :: "real ⇒ real"
assumes "m ≤ n"
and der: "⋀x. x ∈ {of_nat m - 1..of_nat n} ⟹ (g has_field_derivative f x) (at x)"
and le:  "⋀x y. ⟦real m - 1 ≤ x; x ≤ y; y ≤ real n⟧ ⟹ f y ≤ f x"
shows "(∑k = m..n. f (of_nat k)) ≤ g (of_nat n) - g (of_nat m - 1)"
proof -
have "(∑k = m..n. f (of_nat k)) ≤ (∑k = m..n. g (of_nat(Suc k) - 1) - g (of_nat k - 1))"
proof (rule sum_mono, clarsimp)
fix r
assume r: "m ≤ r" "r ≤ n"
hence "∃z>real r - 1. z < real r ∧ g (real r) - g (real r - 1) = (real r - (real r - 1)) * f z"
using assms by (intro MVT2) auto
hence "∃z∈{of_nat r - 1..of_nat r}. g (real r) - g (real r - 1) = f z" by auto
then obtain u::real where u: "u ∈ {of_nat r - 1..of_nat r}"
and eq: "g r - g (of_nat r - 1) = f u" by blast
have "real m ≤ u + 1"
using r u by auto
then have "f (of_nat r) ≤ f u"
using r(2) and u by (intro le) auto
then show "f (of_nat r) ≤ g r - g (of_nat r - 1)"
qed
also have "… ≤ g (of_nat n) - g (of_nat m - 1)"
using ‹m ≤ n› by (subst sum_Suc_diff) auto
finally show ?thesis .
qed

lemma Mertens_lemma:
assumes "n ≠ 0"
shows "¦(∑d = 1..n. mangoldt d / real d) - ln n¦ ≤ 4"
proof -
have *: "⟦abs(s' - nl + n) ≤ a; abs(s' - s) ≤ (k - 1) * n - a⟧
⟹ abs(s - nl) ≤ n * k" for s' s k nl a::real
by (auto simp: algebra_simps abs_if split: if_split_asm)
have le: "¦(∑d=1..n. mangoldt d * floor (n / d)) - n * ln n + n¦ ≤ 1 + ln n"
using ln_fact_bounds ln_fact_conv_mangoldt assms by simp
have "¦real n * ((∑d = 1..n. mangoldt d / real d) - ln n)¦ =
¦((∑d = 1..n. real n * mangoldt d / real d) - n * ln n)¦"
also have "… ≤ real n * 4"
proof (rule * [OF le])
have "¦(∑d = 1..n. mangoldt d * ⌊n / d⌋) - (∑d = 1..n. n * mangoldt d / d)¦
= ¦∑d = 1..n. mangoldt d * (⌊n / d⌋ - n / d)¦"
also have "… ≤ psi n" (is "¦?sm¦ ≤ ?rhs")
proof -
have "-?sm = (∑d = 1..n. mangoldt d * (n/d - ⌊n/d⌋))"
also have "… ≤ (∑d = 1..n. mangoldt d * 1)"
by (intro sum_mono mult_left_mono mangoldt_nonneg) linarith+
finally have "-?sm ≤ ?rhs" by (simp add: psi_def)
moreover
have "?sm ≤ 0"
using mangoldt_nonneg by (simp add: mult_le_0_iff sum_nonpos)
ultimately show ?thesis by (simp add: abs_if)
qed
also have "… ≤ 3/2 * real n"
by (rule psi_ubound_3_2)
also have "…≤ (4 - 1) * real n - (1 + ln n)"
using ln_le_minus_one [of n] assms by (simp add: divide_simps)
finally
show "¦(∑d = 1..n. mangoldt d * real_of_int ⌊real n / real d⌋) -
(∑d = 1..n. real n * mangoldt d / real d)¦
≤ (4 - 1) * real n - (1 + ln n)" .
qed
finally have "¦real n * ((∑d = 1..n. mangoldt d / real d) - ln n)¦ ≤ real n * 4" .
then show ?thesis
using assms mult_le_cancel_left_pos by (simp add: abs_mult)
qed

lemma Mertens_mangoldt_versus_ln:
assumes "I ⊆ {1..n}"
shows "¦(∑i∈I. mangoldt i / i) - (∑p | prime p ∧ p ∈ I. ln p / p)¦ ≤ 3"
(is "¦?lhs¦ ≤ 3")
proof (cases "n = 0")
case True
with assms show ?thesis by simp
next
case False
have "finite I"
using assms finite_subset by blast
have "0 ≤ (∑i∈I. mangoldt i / i - (if prime i then ln i / i else 0))"
using mangoldt_nonneg by (intro sum_nonneg) simp_all
moreover have "… ≤ (∑i = 1..n. mangoldt i / i - (if prime i then ln i / i else 0))"
using assms by (intro sum_mono2) (auto simp: mangoldt_nonneg)
ultimately have *: "¦∑i∈I. mangoldt i / i - (if prime i then ln i / i else 0)¦
≤ ¦∑i = 1..n. mangoldt i / i - (if prime i then ln i / i else 0)¦"
by linarith
moreover have "?lhs = (∑i∈I. mangoldt i / i - (if prime i then ln i / i else 0))"
"(∑i = 1..n. mangoldt i / i - (if prime i then ln i / i else 0))
= (∑d = 1..n. mangoldt d / d) - (∑p | prime p ∧ p ∈ {1..n}. ln p / p)"
using sum.inter_restrict [of _ "λi. ln (real i) / i" "Collect prime", symmetric]
by (force simp: sum_subtractf ‹finite I› intro: sum.cong)+
ultimately have "¦?lhs¦ ≤ ¦(∑d = 1..n. mangoldt d / d) -
(∑p | prime p ∧ p ∈ {1..n}. ln p / p)¦" by linarith
also have "… ≤ 3"
proof -
have eq_sm: "(∑i = 1..n. mangoldt i / i) =
(∑i ∈ {p^k |p k. prime p ∧ p^k ≤ n ∧ k ≥ 1}. mangoldt i / i)"
proof (intro sum.mono_neutral_right ballI, goal_cases)
case (3 i)
hence "¬primepow i" by (auto simp: primepow_def Suc_le_eq)
thus ?case by (simp add: mangoldt_def)
qed (auto simp: Suc_le_eq prime_gt_0_nat)
have "(∑i = 1..n. mangoldt i / i) - (∑p | prime p ∧ p ∈ {1..n}. ln p / p) =
(∑i ∈ {p^k |p k. prime p ∧ p^k ≤ n ∧ k ≥ 2}. mangoldt i / i)"
proof -
have eq: "{p ^ k |p k. prime p ∧ p ^ k ≤ n ∧ 1 ≤ k} =
{p ^ k |p k. prime p ∧ p ^ k ≤ n ∧ 2 ≤ k} ∪ {p. prime p ∧ p ∈ {1..n}}"
(is "?A = ?B ∪ ?C")
proof (intro equalityI subsetI; (elim UnE)?)
fix x assume "x ∈ ?A"
then obtain p k where "x = p ^ k" "prime p" "p ^ k ≤ n" "k ≥ 1" by auto
thus "x ∈ ?B ∪ ?C"
by (cases "k ≥ 2") (auto simp: prime_power_iff Suc_le_eq)
next
fix x assume "x ∈ ?B"
then obtain p k where "x = p ^ k" "prime p" "p ^ k ≤ n" "k ≥ 1" by auto
thus "x ∈ ?A" by (auto simp: prime_power_iff Suc_le_eq)
next
fix x assume "x ∈ ?C"
then obtain p where "x = p ^ 1" "1 ≥ (1::nat)" "prime p" "p ^ 1 ≤ n" by auto
thus "x ∈ ?A" by blast
qed
have eqln: "(∑p | prime p ∧ p ∈ {1..n}. ln p / p) =
(∑p | prime p ∧ p ∈ {1..n}. mangoldt p / p)"
by (rule sum.cong) auto
have "(∑i ∈ {p^k |p k. prime p ∧ p^k ≤ n ∧ k ≥ 1}. mangoldt i / i) =
(∑i ∈ {p ^ k |p k. prime p ∧ p ^ k ≤ n ∧ 2 ≤ k} ∪
{p. prime p ∧ p ∈ {1..n}}. mangoldt i / i)" by (subst eq) simp_all
also have "… = (∑i ∈ {p^k |p k. prime p ∧ p^k ≤ n ∧ k ≥ 2}. mangoldt i / i)
+ (∑p | prime p ∧ p ∈ {1..n}. mangoldt p / p)"
by (intro sum.union_disjoint) (auto simp: prime_power_iff finite_nat_set_iff_bounded_le)
also have "… = (∑i ∈ {p^k |p k. prime p ∧ p^k ≤ n ∧ k ≥ 2}. mangoldt i / i)
+ (∑p | prime p ∧ p ∈ {1..n}. ln p / p)" by (simp only: eqln)
finally show ?thesis
using eq_sm by auto
qed
have "(∑p | prime p ∧ p ∈ {1..n}. ln p / p) ≤ (∑p | prime p ∧ p ∈ {1..n}. mangoldt p / p)"
using mangoldt_nonneg by (auto intro: sum_mono)
also have "… ≤ (∑i = Suc 0..n. mangoldt i / i)"
by (intro sum_mono2) (auto simp: mangoldt_nonneg)
finally have "0 ≤ (∑i = 1..n. mangoldt i / i) - (∑p | prime p ∧ p ∈ {1..n}. ln p / p)"
by simp
moreover have "(∑i = 1..n. mangoldt i / i) - (∑p | prime p ∧ p ∈ {1..n}. ln p / p) ≤ 3"
(is "?M - ?L ≤ 3")
proof -
have *: "∃q. ∃j∈{1..n}. prime q ∧ 1 ≤ q ∧ q ≤ n ∧
(q ^ j = p ^ k ∧ mangoldt (p ^ k) / real p ^ k ≤ ln (real q) / real q ^ j)"
if "prime p" "p ^ k ≤ n" "1 ≤ k" for p k
proof -
have "mangoldt (p ^ k) / real p ^ k ≤ ln p / p ^ k"
using that by (simp add: divide_simps)
moreover have "p ≤ n"
using that self_le_power[of p k] by (simp add: prime_ge_Suc_0_nat)
moreover have "k ≤ n"
proof -
have "k < 2^k"
using of_nat_less_two_power of_nat_less_numeral_power_cancel_iff by blast
also have "… ≤ p^k"
by (simp add: power_mono prime_ge_2_nat that)
also have "… ≤ n"
finally show ?thesis by (simp add: that)
qed
ultimately show ?thesis
using prime_ge_1_nat that by auto (use atLeastAtMost_iff in blast)
qed
have finite: "finite {p ^ k |p k. prime p ∧ p ^ k ≤ n ∧ 1 ≤ k}"
by (rule finite_subset[of _ "{..n}"]) auto
have "?M ≤ (∑(x, k)∈{p. prime p ∧ p ∈ {1..n}} × {1..n}. ln (real x) / real x ^ k)"
by (subst eq_sm, intro sum_le_included [where i = "λ(p,k). p^k"])
(insert * finite, auto)
also have "… = (∑p | prime p ∧ p ∈ {1..n}. (∑k = 1..n. ln p / p^k))"
by (subst sum.Sigma) auto
also have "… = ?L + (∑p | prime p ∧ p ∈ {1..n}. (∑k = 2..n. ln p / p^k))"
finally have "?M - ?L ≤ (∑p | prime p ∧ p ∈ {1..n}. (∑k = 2..n. ln p / p^k))"
also have "… = (∑p | prime p ∧ p ∈ {1..n}. ln p * (∑k = 2..n. inverse p ^ k))"
also have "… = (∑p | prime p ∧ p ∈ {1..n}.
ln p * (((inverse p)⇧2 - inverse p ^ Suc n) / (1 - inverse p)))"
by (intro sum.cong refl) (simp add: sum_gp)
also have "… ≤ (∑p | prime p ∧ p ∈ {1..n}. ln p * inverse (real (p * (p - 1))))"
by (intro sum_mono mult_left_mono)
(auto simp: divide_simps power2_eq_square of_nat_diff mult_less_0_iff)
also have "… ≤ (∑p = 2..n. ln p * inverse (real (p * (p - 1))))"
by (rule sum_mono2) (use prime_ge_2_nat in auto)
also have "… ≤ (∑i = 2..n. ln i / (i - 1)⇧2)"
unfolding divide_inverse power2_eq_square mult.assoc
by (auto intro: sum_mono mult_left_mono mult_right_mono)
also have "… ≤ 3"
proof (cases "n ≥ 3")
case False then show ?thesis
proof (cases "n ≥ 2")
case False then show ?thesis by simp
next
case True
then have "n = 2" using False by linarith
with ln_le_minus_one [of 2] show ?thesis by simp
qed
next
case True
have "(∑i = 3..n. ln (real i) / (real (i - Suc 0))⇧2)
≤ (ln (of_nat n - 1)) - (ln (of_nat n)) - (ln (of_nat n) / (of_nat n - 1)) + 2 * ln 2"
proof -
have 1: "((λz. ln (z - 1) - ln z - ln z / (z - 1)) has_field_derivative ln x / (x - 1)⇧2) (at x)"
if x: "x ∈ {2..real n}" for x
by (rule derivative_eq_intros | rule refl |
(use x in ‹force simp: power2_eq_square divide_simps›))+
have 2: "ln y / (y - 1)⇧2 ≤ ln x / (x - 1)⇧2" if xy: "2 ≤ x" "x ≤ y" "y ≤ real n" for x y
proof (cases "x = y")
case False
define f' :: "real ⇒ real"
where "f' = (λu. ((u - 1)⇧2 / u - ln u * (2 * u - 2)) / (u - 1) ^ 4)"
have f'_altdef: "f' u = inverse u * inverse ((u - 1)⇧2) - 2 * ln u / (u - 1) ^ 3"
if u: "u ∈ {x..y}" for u::real unfolding f'_def using u (* TODO ugly *)
have deriv: "((λz. ln z / (z - 1)⇧2) has_field_derivative f' u) (at u)"
if u: "u ∈ {x..y}" for u::real unfolding f'_def
by (rule derivative_eq_intros refl | (use u xy in ‹force simp: divide_simps›))+
hence "∃z>x. z < y ∧ ln y / (y - 1)⇧2 - ln x / (x - 1)⇧2 = (y - x) * f' z"
using xy and ‹x ≠ y› by (intro MVT2) auto
then obtain ξ::real where "x < ξ" "ξ < y"
and ξ: "ln y / (y - 1)⇧2 - ln x / (x - 1)⇧2 = (y - x) * f' ξ" by blast
have "f' ξ ≤ 0"
proof -
have "2/3 ≤ ln (2::real)" by (fact ln_2_ge')
also have "… ≤ ln ξ"
using ‹x < ξ› xy by auto
finally have "1 ≤ 2 * ln ξ" by simp
then have *: "ξ ≤ ξ * (2 * ln ξ)"
using ‹x < ξ› xy by auto
hence "ξ - 1 ≤ ln ξ * 2 * ξ" by (simp add: algebra_simps)
hence "1 / (ξ * (ξ - 1)⇧2) ≤ ln ξ * 2 / (ξ - 1) ^ 3"
using xy ‹x < ξ› by (simp add: divide_simps power_eq_if)
thus ?thesis using xy ‹x < ξ› ‹ξ < y› by (subst f'_altdef) (auto simp: divide_simps)
qed
then have "(ln y / (y - 1)⇧2 - ln x / (x - 1)⇧2) ≤ 0"
using ‹x ≤ y› by (simp add: mult_le_0_iff ξ)
then show ?thesis by simp
qed simp_all
show ?thesis
using sum_integral_ubound_decreasing'
[OF ‹3 ≤ n›, of "λz. ln(z-1) - ln z - ln z / (z - 1)" "λz. ln z / (z-1)⇧2"]
1 2 ‹3 ≤ n›
by (auto simp: in_Reals_norm of_nat_diff)
qed
also have "… ≤ 2"
proof -
have "ln (real n - 1) - ln n ≤ 0"  "0 ≤ ln n / (real n - 1)"
using ‹3 ≤ n› by auto
then have "ln (real n - 1) - ln n - ln n / (real n - 1) ≤ 0"
by linarith
with ln_2_less_1 show ?thesis by linarith
qed
also have "… ≤ 3 - ln 2"
using ln_2_less_1 by (simp add: algebra_simps)
finally show ?thesis
using True by (simp add: algebra_simps sum.atLeast_Suc_atMost [of 2 n])
qed
finally show ?thesis .
qed
ultimately show ?thesis
by linarith
qed
finally show ?thesis .
qed

proposition Mertens:
assumes "n ≠ 0"
shows "¦(∑p | prime p ∧ p ≤ n. ln p / of_nat p) - ln n¦ ≤ 7"
proof -
have "¦(∑d = 1..n. mangoldt d / real d) - (∑p | prime p ∧ p ∈ {1..n}. ln (real p) / real p)¦
≤ 7 - 4" using Mertens_mangoldt_versus_ln [of "{1..n}" n] by simp_all
also have "{p. prime p ∧ p ∈ {1..n}} = {p. prime p ∧ p ≤ n}"
using atLeastAtMost_iff prime_ge_1_nat by blast
finally have "¦(∑d = 1..n. mangoldt d / real d) - (∑p∈…. ln (real p) / real p)¦ ≤ 7 - 4" .
moreover from assms have "¦(∑d = 1..n. mangoldt d / real d) - ln n¦ ≤ 4"
by (rule Mertens_lemma)
ultimately show ?thesis by linarith
qed

end


# File ‹bertrand.ML›

signature BERTRAND =
sig

type cache = (int * thm) list
datatype primepow_thm = PrimepowThm of thm * thm | NotPrimepowThm of thm

val prime_conv : Proof.context -> cache -> conv
val primepow_conv : Proof.context -> cache -> conv
val pre_mangoldt_conv : Proof.context -> cache -> conv

val prove_primepow : Proof.context -> cache -> term -> primepow_thm
val prove_pre_mangoldt : Proof.context -> cache -> term -> thm
val prove_psi : Proof.context -> int -> (int * int * thm) list

val mk_prime_cache : Proof.context -> int -> cache

end

structure Bertrand : BERTRAND =
struct

type cache = (int * thm) list

datatype primepow_cert = Primepow of int * int | NotPrimepow of int * int
datatype primepow_thm = PrimepowThm of thm * thm | NotPrimepowThm of thm

val mk_nat = HOLogic.mk_number @{typ nat}

fun mk_primepow_cert n =
let
fun find_prime_divisor n =
let fun go k = if n mod k = 0 then k else go (k + 1) in go 2 end
fun divide_out p m acc =
if m mod p = 0 then divide_out p (m div p) (acc + 1) else (acc, m)
val p = find_prime_divisor n
val (k, m) = divide_out p n 0
in
if m = 1 then Primepow (p, k) else NotPrimepow (p, find_prime_divisor m)
end

fun bool_thm_to_eq_thm thm =
case HOLogic.dest_Trueprop (Thm.concl_of thm) of
@{term "HOL.Not"} $_ => thm RS @{thm Eq_FalseI} | _ => thm RS @{thm Eq_TrueI} fun prime_conv ctxt cache ct = case Thm.term_of ct of @{const "HOL.Trueprop"}$ (@{term "prime :: nat ⇒ bool"} $p) => Pratt.prove_prime cache (snd (HOLogic.dest_number p)) ctxt |> fst |> the |> bool_thm_to_eq_thm | _ => raise CTERM ("prime_conv", [ct]) fun prime_sieve n = let fun go ps k = if k > n then rev ps else if exists (fn p => k mod p = 0) ps then go ps (k + 1) else go (k :: ps) (k + 1) in go [] 2 end fun fold_accum f xs acc = fold (fn x => fn (ys, acc) => case f x acc of (y, acc') => (y :: ys, acc')) xs ([], acc) |> apfst rev fun mk_prime_cache ctxt n = let val ps = prime_sieve n in fold_accum (fn p => fn cache => Pratt.prove_prime cache p ctxt) ps [] |> snd end fun prove_primepow ctxt cache t = let val (_, n) = HOLogic.dest_number t fun inst xs = Drule.infer_instantiate' ctxt (map (SOME o Thm.cterm_of ctxt o mk_nat) xs) val prove = Goal.prove ctxt [] [] fun prove_prime' p = the (fst (Pratt.prove_prime cache p ctxt)) in if n <= 0 then raise TERM ("prove_primepow", [t]) else if n = 1 then NotPrimepowThm @{thm not_primepow_1_nat} else case mk_primepow_cert n of Primepow (p, k) => let val thm = prove_prime' p RS inst [p, k, n] @{thm primepowI} val thm' = prove (Thm.concl_of thm) (fn {context = ctxt, ...} => HEADGOAL (resolve_tac ctxt [thm]) THEN ALLGOALS (Simplifier.simp_tac ctxt)) in PrimepowThm (thm' RS @{thm conjunct1}, thm' RS @{thm conjunct2}) end | NotPrimepow (p, q) => let val thm = inst [p, q, n] @{thm not_primepowI} OF (map prove_prime' [p, q]) val thm' = prove (Thm.concl_of thm) (fn {context = ctxt, ...} => HEADGOAL (resolve_tac ctxt [thm]) THEN ALLGOALS (Simplifier.simp_tac ctxt)) in NotPrimepowThm thm' end end fun primepow_conv ctxt cache ct = case prove_primepow ctxt cache (Thm.term_of ct) of PrimepowThm (thm, _) => bool_thm_to_eq_thm thm | NotPrimepowThm thm => bool_thm_to_eq_thm thm fun prove_pre_mangoldt ctxt cache t = case prove_primepow ctxt cache t of PrimepowThm (thm1, thm2) => @{thm pre_mangoldt_primepow} OF [thm1, thm2] | NotPrimepowThm thm => thm RS @{thm pre_mangoldt_notprimepow} fun pre_mangoldt_conv ctxt cache = (fn thm => thm RS @{thm eq_reflection}) o prove_pre_mangoldt ctxt cache o Thm.term_of val nat_eq_ss = simpset_of (put_simpset HOL_basic_ss @{context} addsimps @{thms Numeral_Simprocs.semiring_norm}) fun prove_nat_eq ctxt (t1, t2) = let val goal = HOLogic.mk_eq (t1, t2) |> HOLogic.mk_Trueprop fun tac {context = ctxt, ...} = HEADGOAL (resolve_tac ctxt @{thms mult_1_left mult_1_right}) ORELSE HEADGOAL (Simplifier.simp_tac (put_simpset nat_eq_ss ctxt)) in Goal.prove ctxt [] [] goal tac end fun prove_psi ctxt n = let val cache = mk_prime_cache ctxt n fun go thm x x' k acc = if k > n then rev acc else let val pre_mangoldt_thm = prove_pre_mangoldt ctxt cache (mk_nat k) val y = pre_mangoldt_thm |> Thm.concl_of |> HOLogic.dest_Trueprop |> HOLogic.dest_eq |> snd val y' = HOLogic.dest_number y |> snd val eq_thm1 = prove_nat_eq ctxt (@{term "(+) :: nat ⇒ _"}$ mk_nat (k - 1) $@{term "1 :: nat"}, mk_nat k) val z = if y' = 1 then x else mk_nat (x' * y') val eq_thm2 = prove_nat_eq ctxt (@{term "(*) :: nat => _"}$ x \$ y, z)
val thm' = @{thm eval_psi_aux2} OF [thm, pre_mangoldt_thm, eq_thm1, eq_thm2]
in
go thm' z (x' * y') (k + 1) ((k - 1, x', thm) :: acc)
end
in
go @{thm eval_psi_aux1} @{term "numeral Num.One :: nat"} 1 1 []
end

end
`