Session Bertrands_Postulate

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 
        [symmetric] of_nat_add [symmetric] del: of_nat_mult of_nat_add)

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"
    by (simp_all add: aprimedivisor_primepow prime_elem_multiplicity_mult_distrib)
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"
    by (simp_all add: aprimedivisor_primepow prime_elem_multiplicity_mult_distrib)
  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)"
    by (simp_all add: aprimedivisor_primepow prime_elem_multiplicity_mult_distrib)
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"
    by (simp_all add: aprimedivisor_primepow prime_elem_multiplicity_mult_distrib)
  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"
      by (simp add: sum.distrib)
    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))"
      by (simp add: sum.distrib)
    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)"
      by (simp add: algebra_simps)
    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"
      by (simp add: algebra_simps mult_left_mono)
    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))"
        by (intro ln_add_one_self_le_self) simp_all
      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)"
      by (simp add: ln_div)
    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))"
      by (simp add: algebra_simps)
    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"
    by (simp add: algebra_simps)
  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"
    by (simp add: algebra_simps)
  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 *: "(dprimepow_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)"
      by (simp add: of_nat_prod_mset)
  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"
        by (simp add: floor_divide_of_nat_eq)
      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)"
    by (simp add: sum_distrib_left mult_ac)
  also have "(d = 1..n. mangoldt d * real_of_int real n / real d) -  = 
               ln (fact n) - 2 * ln (fact (n div 2))"
    by (simp add: ln_fact_conv_mangoldt)
  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 "xinsert y A. P x"
  shows   "P y" "xA. 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))"
  by (simp add: psi_def)

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"
    by (simp add: field_simps log_def)
  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))" 
    by (simp add: powr_realpow [symmetric])
  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" 
      by (simp add: algebra_simps)
    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"
    by (simp add: field_simps)
  finally show ?thesis using c  0 by (simp_all add: mult_left_mono)
qed

end  

lemma psi_ubound_log_double_cases:
  assumes "nm. 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   "nm'. 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:
  "n1024. psi n  551 / 256 * ln 2 * real n"
proof -
  from psi_ubound_log_128 have "n128. psi n  3 / 2 * ln 2 * real n" by simp
  hence "n256. 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 "n512. 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 "n1024. 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"
        by (simp add: left_diff_distrib)
      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"
        by (simp add: add_pos_pos)
      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)
  qed (simp_all add: of_nat_ge_1_iff Suc_le_eq)
  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))"
          by (simp add: aprimedivisor_primepow_power)
        also from True have "ln () = mangoldt (Suc (Discrete.sqrt n))"
          by (simp add: mangoldt_def)
        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))"
          by (simp add: mangoldt_def)
        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"] 
  by (simp add: psi_odd_def)

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"
    by (simp add: minus_mod_eq_div_mult[symmetric])
  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))"
    by (simp add: psi_def)
  also have " = (d = 1..n. e = 1..n div d. (-1)^(d+1) * mangoldt e)"
    by (simp add: sum_distrib_left)
  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)"
    by (simp add: mult_ac div_mult2_eq)
  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)"
    by (simp add:  ln_fact_conv_mangoldt)
  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)"
    by (simp add: sum_distrib_left mult_ac)
  also have "(d=1..n. mangoldt d * n / d) -  = 
               (d=1..n. (mangoldt d * n / d - mangoldt d * 2 * real (n div 2) / d))"
    by (simp add: sum_subtractf)
  also have " = (d=1..n. mangoldt d * (n / d - 2 * real (n div 2) / d))"
    by (simp add: algebra_simps)
  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)"
    by (simp add: mult_ac)
  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))"
    by (simp add: psi_mono div_le_mono2)
  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))"
    by (simp add: psi_mono div_le_mono2)
  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)"
    by (simp add: sum.atLeast_Suc_atMost numeral_2_eq_2)
  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"
    by (simp add: power2_eq_square)
  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 "pprime_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)"
      by (simp add: eq)
  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)¦"
    by (simp add: algebra_simps sum_distrib_left)
  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)¦"
      by (simp add: sum_subtractf algebra_simps)
    also have "  psi n" (is "¦?sm¦  ?rhs")
    proof -
      have "-?sm = (d = 1..n. mangoldt d * (n/d - n/d))"
        by (simp add: sum_subtractf algebra_simps)
      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 "¦(iI. 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  (iI. 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 *: "¦iI. 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 = (iI. 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"
              by (simp add: that)
            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))"
          by (simp add: comm_monoid_add_class.sum.distrib sum.atLeast_Suc_atMost numeral_2_eq_2)
        finally have "?M - ?L  (p | prime p  p  {1..n}. (k = 2..n. ln p / p^k))"
          by (simp add: algebra_simps)
        also have " = (p | prime p  p  {1..n}. ln p * (k = 2..n. inverse p ^ k))"
          by (simp add: field_simps sum_distrib_left)
        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 *)
                by (simp add: eval_nat_numeral divide_simps) (simp add: algebra_simps)?
              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