Theory Bertrand

(*
  File:     Bertrand.thy
  Authors:  Julian Biendarra, Manuel Eberl <manuel@pruvisto.org>, 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_Functions"
    "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_pos 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 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 lthy =>
let
  fun tac ctxt =
    let
      val psi_cache = Bertrand.prove_psi ctxt 129
      fun prove_psi_ineqs ctxt =
        let
          fun tac goal_ctxt = 
            HEADGOAL (resolve_tac goal_ctxt @{thms eval_psi_ineq_aux2} THEN'
              Simplifier.simp_tac goal_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 o #context)
            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 (fn {context = goal_ctxt, ...} =>
                  HEADGOAL (Simplifier.simp_tac goal_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
                  val thm' =
                    Goal.prove ctxt [] [] (Thm.cprem_of thm 1 |> Thm.term_of)
                      (fn {context = goal_ctxt, ...} =>
                        HEADGOAL (Simplifier.simp_tac goal_ctxt))
                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 lthy [] [] @{prop "n{0..128}. psi n  3 / 2 * ln 2 * n"} (tac o #context)
in
  Local_Theory.note ((@{binding psi_ubound_log_128}, []), [thm]) lthy |> 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:<