Theory Chebyshev_Prime_Bounds

theory Chebyshev_Prime_Bounds
imports
  "Prime_Number_Theorem.Prime_Counting_Functions"
  "Prime_Distribution_Elementary.Prime_Distribution_Elementary_Library"
  "Prime_Distribution_Elementary.Primorial"
  "HOL-Decision_Procs.Approximation"
  "HOL-Library.Code_Target_Numeral"
  Chebyshev_Prime_Exhaust
begin

subsection ‹Auxiliary material›

context comm_monoid_set
begin

lemma union_disjoint':
  assumes "finite C" "A  B = C" "A  B = {}"
  shows   "f (F g A) (F g B) = F g C"
  using union_disjoint[of A B g] assms by auto

end

lemma sum_mset_nonneg:
  fixes X :: "'a :: ordered_comm_monoid_add multiset"
  shows "(x. x ∈# X  x  0)  sum_mset X  0"
  by (induction X) (auto)

lemma of_int_sum_mset: "of_int (sum_mset M) = sum_mset (image_mset of_int M)"
  by (induction M) auto

lemma sum_sum_mset: "(xA. y∈#B. f x y) = (y∈#B. xA. f x y)"
  by (induction B) (auto simp: algebra_simps sum.distrib)

lemma sum_mset_diff_distrib:
  fixes f g :: "'a  'b :: ab_group_add"
  shows   "(x∈#A. f x - g x) = (x∈#A. f x) - (x∈#A. g x)"
  by (induction A) (auto simp: algebra_simps)

lemma sum_mset_neg_distrib:
  fixes f :: "'a  'b :: ab_group_add"
  shows   "(x∈#A. -f x) = -(x∈#A. f x)"
  by (induction A) (auto simp: algebra_simps)


subsection ‹Bounds for the remainder in Stirling's approximation›

definition ln_fact_remainder :: "real  real" where
  "ln_fact_remainder x = ln (fact (nat x)) - (x * ln x - x)"

lemma ln_fact_remainder_bounds:
  assumes x: "x  3"
  shows "ln_fact_remainder x   ln x / 2 + ln (2 * pi) / 2 + 1 / (12 * x)"
    and "ln_fact_remainder x  -ln x / 2 + ln (2 * pi) / 2 - 1 / (2 * x)"
proof -
  define n where "n = nat x"
  define f where "f = (λt. t * (ln t - 1) + ln t / 2 :: real)"

  have "ln x  1"
  proof -
    have "1  ln (3 :: real)"
      by (approximation 10)
    also have "ln 3  ln x"
      using assms by simp
    finally show ?thesis .
  qed 

  have n: "n  1"
    using x by (auto simp: n_def le_nat_iff)
  have "ln_fact_remainder x = ln (fact n) + ln x / 2 - f x"
    by (simp add: ln_fact_remainder_def n_def f_def algebra_simps)
  also have "ln (fact n)  ln (2 * pi) / 2 + f n + 1 / (12 * n)"
    using ln_fact_bounds(2)[of n] n by (auto simp: f_def ln_mult add_divide_distrib algebra_simps)
  also have " + ln x / 2 - f x = ln x / 2 + (f n - f x) + ln (2 * pi) / 2 + 1 / (12 * n)"
    using n by (simp add: algebra_simps ln_mult)
  also have "f n  f x"
    unfolding f_def using assms ln x  1
    by (intro add_mono mult_mono) (auto simp: n_def)
  finally show "ln_fact_remainder x  ln x / 2 + ln (2 * pi) / 2 + 1 / (12 * x)"
    using assms by (simp add: n_def)

  define f' :: "real  real" where "f' = (λx. ln x + 1 / (2 * x))"
  have f'_mono: "f' x  f' y" if "x  y" "x  1 / 2" for x y :: real
    using that(1)
  proof (rule DERIV_nonneg_imp_nondecreasing)
    fix t assume t: "t  x" "t  y"
    hence "t > 0"
      using x  1 / 2 by auto
    have "(t - 1 / 2) / t ^ 2  0"
      using t that by auto
    have "(f' has_field_derivative (1 / t - 1 / (2 * t ^ 2))) (at t)"
      using t > 0 by (auto simp: f'_def power2_eq_square intro!: derivative_eq_intros)
    also have "1 / t - 1 / (2 * t ^ 2) = (t - 1 / 2) / t ^ 2"
      using t > 0 by (simp add: field_simps eval_nat_numeral del: div_diff)
    finally show "y. (f' has_real_derivative y) (at t)  0  y"
      using (t - 1 / 2) / t ^ 2  0 by blast
  qed

  have f'_nonneg: "f' t  0" if "t  3" for t
  proof -
    have "0  f' 3"
      unfolding f'_def by (approximation 10)
    also have "f' 3  f' t"
      by (rule f'_mono) (use that in auto)
    finally show ?thesis .
  qed

  have "f x - f n  f' x * frac x"
  proof (cases "n < x")
    case False
    hence "x = n"
      using assms unfolding n_def by linarith
    thus ?thesis using f'_nonneg[of x] assms
      by (simp add: n_def)
  next
    case True
    have "z::real. z > n  z < x  f x - f n = (x - n) * f' z"
      using True assms n
      by (intro MVT2) (auto intro!: derivative_eq_intros simp: f_def f'_def)
    then obtain z :: real where z: "z > n" "z < x" "f x - f n = (x - n) * f' z"
      by blast
    have "f' z  f' x"
      by (rule f'_mono) (use z assms n in auto)

    have "f x - f n = (x - n) * f' z"
      by fact
    also have "  (x - n) * f' x"
      using f' z  f' x True by (intro mult_left_mono ) auto
    also have "x - n = frac x"
      using assms by (simp add: n_def frac_def)
    finally show ?thesis
      by (simp add: mult_ac)
  qed

  also have "  f' x * 1"
    using frac_lt_1[of x] f'_nonneg[of x] assms
    by (intro mult_left_mono) auto
  finally have "f n - f x  -1 / (2 * x) - ln x"
    by (simp add: f'_def)

  have "-ln x / 2 - 1 / (2 * x) + ln (2 * pi) / 2 = 
        ln x / 2 + (-1 / (2 * x) - ln x) + ln (2 * pi) / 2"
    by (simp add: algebra_simps)
  also have "-1 / (2 * x) - ln x  f n - f x"
    by fact
  also have "ln x / 2 + (f n - f x) + ln (2 * pi) / 2 =
          ln (2 * pi) / 2 + f n + ln x / 2 - f x"
    by (simp add: algebra_simps)
  also have "ln (2 * pi) / 2 + f n  ln (fact n)"
    using ln_fact_bounds(1)[of n] n by (auto simp: f_def ln_mult add_divide_distrib algebra_simps)
  also have "ln (fact n) + ln x / 2 - f x = ln_fact_remainder x"
    by (simp add: ln_fact_remainder_def f_def n_def algebra_simps)
  finally show "ln_fact_remainder x  -ln x / 2 + ln (2 * pi) / 2 - 1 / (2 * x)"
    by simp
qed

lemma abs_ln_fact_remainder_bounds:
  assumes x: "x  3"
  shows "¦ln_fact_remainder x¦ < ln x / 2 + 1"
proof -
  have "ln_fact_remainder x  ln x / 2 + (ln (2 * pi) / 2 + 1 / (12 * x))"
    using ln_fact_remainder_bounds(1)[of x] assms by (simp add: algebra_simps)
  also have "1 / (12 * x)  1 / 36"
    using assms by auto
  also have "ln (2 * pi) / 2 + 1 / 36 < 1"
    by (approximation 10)
  finally have less: "ln_fact_remainder x < ln x / 2 + 1"
    by simp

  have "-(ln x / 2 + 1) = -ln x / 2 + (-1)"
    by simp
  also have "-1 < 0 - 1 / (2 * x)"
    using assms by simp
  also have "0  ln (2 * pi) / 2"
    using pi_gt3 by simp
  also have "-ln x / 2 + (ln (2 * pi) / 2 - 1 / (2 * x))  ln_fact_remainder x"
    using ln_fact_remainder_bounds(2)[of x] assms by (simp add: algebra_simps)
  finally have "- (ln x / 2 + 1) < ln_fact_remainder x" by - simp_all
  with less show ?thesis
    by linarith
qed


subsection ‹Approximating ψ›

unbundle prime_counting_syntax

lemma primes_psi_lower_rec:
  fixes f :: "real  real"
  assumes "x. x  x0  f x  f (x / c) + h x"
  assumes  "x0 > 0" "x * c  x0 * c ^ n" "c  1"
  shows   "f x  f (x / c ^ n) + (k<n. h (x / c ^ k))"
  using assms(2-)
proof (induction n arbitrary: x)
  case 0
  thus ?case by auto
next
  case (Suc n)
  have "0 < x0 * c ^ n"
    using Suc.prems by auto
  also have "  x"
    using Suc.prems by auto
  finally have "x > 0" .

  have "x0 * c ^ n  1 * x"
    using Suc.prems by simp
  also have "1 * x  c * x"
    by (rule mult_right_mono) (use Suc.prems x > 0 in auto)
  finally have "f x  f (x / c ^ n) + (k<n. h (x / c ^ k))"
    by (intro Suc.IH) (use Suc.prems in auto simp: mult_ac)
  also have "f (x / c ^ n)  f (x / c ^ n / c) + h (x / c ^ n)"
    by (rule assms(1)) (use Suc.prems x > 0 in auto simp: field_simps less_imp_le)
  finally show ?case
    by (simp add: mult_ac add_ac)
qed


locale chebyshev_multiset =
  fixes L :: "int multiset"
  assumes L_nonzero: "0 ∉# L"
begin

definition chi_L :: "real  int" ("χL")
  where "chi_L t = (l∈#L. sgn l * t / ¦l¦)"

definition psi_L :: "real  real" ("ψL")
  where "psi_L x = sum_upto (λd. mangoldt d * chi_L (x / d)) x"

definition alpha_L :: real ("αL")
  where "alpha_L = -(l∈#L. ln ¦l¦ / l)"

definition period :: nat
  where "period = nat (Lcm (set_mset L))"

lemma period_pos: "period > 0"
proof -
  have "Lcm (set_mset L)  0"
    using L_nonzero unfolding period_def by (subst Lcm_0_iff) auto
  moreover have "Lcm (set_mset L)  0"
    by auto
  ultimately have "Lcm (set_mset L) > 0"
    by linarith
  thus ?thesis
    by (simp add: period_def)
qed

lemma dvd_period: "l ∈# L  l dvd period"
  unfolding period_def by auto

lemma chi_L_decompose:
  "χL (x + of_int (m * int period)) = χL x + m * int period * (l∈#L. 1 / l)"
proof -
  have "real_of_int (χL (x + of_int (m * int period))) = 
          (l∈#L. of_int (sgn l * (x + of_int m * real period) / real_of_int ¦l¦))"
    by (simp add: chi_L_def of_int_sum_mset multiset.map_comp o_def)
  also have " = (l∈#L. real_of_int (sgn l * (x / of_int ¦l¦)) + m * period / l)"
  proof (intro arg_cong[of _ _ sum_mset] image_mset_cong, goal_cases)
    case (1 l)
    with L_nonzero have [simp]: "l  0"
      by auto
    have "(x + of_int m * real period) / real_of_int ¦l¦ =
           x / of_int ¦l¦ + of_int (m * period div ¦l¦)"
      using dvd_period[of l] 1 by (subst real_of_int_div) (auto simp: field_simps)
    also have "floor  = x / of_int ¦l¦ :: real + m * period div ¦l¦"
      by (subst floor_add_int) auto
    also have "real_of_int  = x / of_int ¦l¦ + m * period / ¦l¦"
      using dvd_period[of l] 1 by (simp add: real_of_int_div)
    also have "sgn l *  = sgn l * x / of_int ¦l¦ + m * period / l"
      by (simp add: sgn_if)
    finally show ?case
      by simp
  qed
  also have " = of_int (χL x) + (l∈#L. m * period / l)"
    by (subst sum_mset.distrib)
       (auto simp: chi_L_def of_int_sum_mset multiset.map_comp o_def)
  also have "(l∈#L. m * period / l) = m * period * (l∈#L. 1 / l)"
    by (simp add: sum_mset_distrib_left)
  finally show ?thesis
    by simp
qed

lemma chi_L_floor: "chi_L (floor x) = chi_L x"
  unfolding chi_L_def
proof (intro arg_cong[of _ _ sum_mset] image_mset_cong, goal_cases)
  case (1 l)
  thus ?case
    using floor_divide_real_eq_div[of "¦l¦" x] floor_divide_of_int_eq[of "x" "¦l¦"]
    by auto
qed  
    
end


locale balanced_chebyshev_multiset = chebyshev_multiset +
  assumes balanced: "(l∈#L. 1 / l) = 0"
begin

lemma chi_L_mod: "χL (of_int (a mod int period)) = χL (of_int a)"
proof -
  have a: "a = a mod period + period * (a div period)"
    by simp
  have "of_int a = real_of_int (a mod int period) +
                   real_of_int (a div int period * int period)"
    by (subst a, unfold of_int_add) auto
  also have "real_of_int (χL ) = real_of_int (χL (real_of_int (a mod int period)))"
    using balanced by (subst chi_L_decompose) auto
  finally show ?thesis
    by linarith
qed

sublocale chi: periodic_fun_simple chi_L "of_int period"
proof
  fix x :: real
  have "χL (x + real_of_int (int period)) = χL (of_int (x + real_of_int (int period) mod int period))"
    unfolding chi_L_mod chi_L_floor ..
  also have "x + real_of_int (int period) mod int period = x mod int period"
    by simp
  also have "χL  = χL x"
    by (simp add: chi_L_mod chi_L_floor)
  finally show "χL (x + real_of_int (int period)) = χL x" .
qed


definition psi_L_remainder where
  "psi_L_remainder x = (l∈#L. sgn l * ln_fact_remainder (x / ¦l¦))"

lemma abs_sum_mset_le:
  fixes f :: "'a  'b :: ordered_ab_group_add_abs"
  shows "¦x∈#A. f x¦  (x∈#A. ¦f x¦)"
  by (induction A) (auto intro: order.trans[OF abs_triangle_ineq])

lemma psi_L_remainder_bounds:
  fixes x :: real
  assumes x: "x  3" "l. l ∈# L  x  3 * ¦l¦"
  shows "¦psi_L_remainder x¦ 
           ln x * size L / 2 - 1/2 * (l∈#L. ln ¦l¦) + size L"
proof -
  have nonzero: "l  0" if "l ∈# L" for l
    using L_nonzero that by auto
  have "psi_L_remainder x = (l∈#L. sgn l * ln_fact_remainder (x / ¦l¦))"
    by (simp add: psi_L_remainder_def)
  also have "¦¦  (l∈#L. ¦sgn l * ln_fact_remainder (x / ¦l¦)¦)"
    by (rule abs_sum_mset_le)
  also have " = (l∈#L. ¦ln_fact_remainder (x / ¦l¦)¦)"
    by (intro arg_cong[of _ _ sum_mset] image_mset_cong)
       (auto simp: nonzero abs_mult simp flip: of_int_abs)
  also have "  (l∈#L. ln (x / ¦l¦) / 2 + 1)"
    using x
    by (intro sum_mset_mono less_imp_le[OF abs_ln_fact_remainder_bounds])
       (auto simp: nonzero field_simps)
  also have " = (l∈#L. 1 / 2 * (ln x - ln ¦l¦) + 1)"
    using assms
    by (intro arg_cong[of _ _ sum_mset] image_mset_cong) (auto simp: algebra_simps ln_div nonzero)
  also have " = ln x / 2 * size L + (-1/2) * (l∈#L. ln ¦l¦) + size L"
    unfolding sum_mset_distrib_left of_int_sum_mset
    by (simp add: sum_mset.distrib sum_mset_diff_distrib diff_divide_distrib sum_mset_neg_distrib)
  finally show ?thesis
    using assms by (simp add: mult_left_mono divide_right_mono add_mono)
qed  

lemma psi_L_eq:
  assumes "x > 0"
  shows   "psi_L x = αL * x + psi_L_remainder x"
proof -
  have "psi_L x = (l∈#L. sgn l *
          sum_upto (λd. mangoldt d * x / (d * ¦l¦)) x)"
    by (simp add: psi_L_def chi_L_def sum_upto_def sum_mset_distrib_left of_int_sum_mset
          multiset.map_comp o_def sum_sum_mset algebra_simps sum_distrib_left sum_distrib_right)
  also have " = (l∈#L. sgn l *
          sum_upto (λd. mangoldt d * x / (d * ¦l¦)) (x / ¦l¦))"
  proof (intro arg_cong[of _ _ sum_mset] image_mset_cong, goal_cases)
    case (1 l)
    have "l  0"
      using 1 L_nonzero by auto

    have "sum_upto (λd. mangoldt d * real_of_int x / real_of_int (int d * ¦l¦)) (x / real_of_int ¦l¦) =
          sum_upto (λd. mangoldt d * real_of_int x / real_of_int (int d * ¦l¦)) x"
      unfolding sum_upto_def
    proof (intro sum.mono_neutral_left subsetI ballI, goal_cases)
      case (2 d)
      hence "real d  x / ¦real_of_int l¦"
        by auto
      also have "  x / 1"
        using l  0 and assms by (intro divide_left_mono) auto
      finally show ?case
        using 2 by auto
    next
      case (3 d)
      hence "x < d * ¦l¦" and "d > 0"
        using l  0 and assms by (auto simp: field_simps)
      hence "x / real_of_int (int d * ¦l¦)  0" and "x / real_of_int (int d * ¦l¦) < 1"
        using assms by auto
      hence "x / real_of_int (int d * ¦l¦) = 0"
        by linarith
      thus ?case
        by simp
    qed auto
    thus ?case
      by simp
  qed

  also have " = (l∈#L. sgn l * ln (fact (nat x/¦l¦)))"
    by (subst ln_fact_conv_sum_mangoldt [symmetric]) (auto simp: mult_ac)

  also have " = (l∈#L. x / l * ln x - x * ln ¦l¦ / l - x / l + sgn l * ln_fact_remainder (x / ¦l¦))"
  proof (intro arg_cong[of _ _ sum_mset] image_mset_cong, goal_cases)
    case (1 l)
    hence [simp]: "l  0"
      using L_nonzero by auto
    have "ln (fact (nat x/¦l¦)) = x / ¦l¦ * ln (x / ¦l¦) - x / ¦l¦ + ln_fact_remainder (x / ¦l¦)"
      by (simp add: ln_fact_remainder_def)
    also have "real_of_int (sgn l) *  = x / l * ln x - x * ln ¦l¦ / l - x / l + sgn l * ln_fact_remainder (x / ¦l¦)"
      using assms by (auto simp: sgn_if ln_div diff_divide_distrib algebra_simps)
    finally show ?case .
  qed
  also have " = (x * ln x - x) * (l∈#L. 1 / l) - x * (l∈#L. ln ¦l¦ / l) + (l∈#L. sgn l * ln_fact_remainder (x / ¦l¦))"
    by (simp add: sum_mset.distrib sum_mset_diff_distrib sum_mset_distrib_left diff_divide_distrib)
  also have " = αL * x + psi_L_remainder x"
    by (subst balanced) (auto simp: alpha_L_def psi_L_remainder_def)
  finally show ?thesis .
qed

lemma primes_psi_lower_bound:
  fixes x C :: real
  defines "x0  Max (insert 3 ((λl. 3 * ¦l¦) ` set_mset L))"
  assumes x: "x  x0"
  assumes chi_le1: "n. n  {0..<period}  χL (real n)  1"
  defines "C  1 / 2 * (l∈#L. ln ¦l¦) - size L"
  shows   "ψ x  αL * x - ln x * size L / 2 + C"
proof -
  have chi_le1': "χL x  1" for x
    proof -
    have "χL x = χL (floor x mod period)"
      by (simp add: chi_L_mod chi_L_floor)
    also have "floor x mod period = real (nat (floor x mod period))"
      using period_pos by auto
    also have "χL   1"
      by (rule chi_le1) (use period_pos in auto simp: nat_less_iff)
    finally show ?thesis .
  qed

  have x0: "x0  3" "l. l ∈# L  x0  3 * ¦l¦"
    unfolding x0_def by auto

  have *: "x * y  x" if "y  1" "x  0" for x y :: real
    using mult_left_mono[OF that] by auto

  have "¦psi_L_remainder x¦  ln x * real (size L) / 2 -
          1 / 2 * (l∈#L. ln (real_of_int ¦l¦)) + real (size L)"
    by (rule psi_L_remainder_bounds)
       (use x x0 in force simp flip: of_int_abs)+
  hence "¦psi_L_remainder x¦  ln x * size L / 2 - C"
    by (simp add: C_def algebra_simps)
  hence "αL * x - ln x * size L / 2 + C  αL * x + psi_L_remainder x"
    by linarith
  also have "αL * x + psi_L_remainder x = ψL x"
    using x x0(1) by (subst psi_L_eq) auto
  also have "ψL x  ψ x"
    unfolding psi_L_def primes_psi_def sum_upto_def
    by (intro sum_mono *) (auto simp: mangoldt_nonneg chi_le1')
  finally show ?thesis
    by (simp add: C_def)
qed

end

lemma psi_lower_bound_precise:
  assumes x: "x  90"
  shows   "ψ x  0.92128 * x - 2.5 * ln x - 1.6"
proof -
  interpret balanced_chebyshev_multiset "{#1, -2, -3, -5, 30#}"
    by unfold_locales auto

  define C :: real where "C = ((ln 2 + (ln 3 + (ln 5 + ln 30))) / 2 - 5)"
  have "alpha_L = ln 2 / 2 - (ln 30 / 30 - ln 5 / 5 - ln 3 / 3)"
    by (simp add: alpha_L_def)
  also have "  0.92128"
    by (approximation 30)
  finally have "alpha_L  0.92128" .
  have "C  -1.6"
    unfolding C_def by (approximation 20)

  have "0.92128 * x - ln x * 5 / 2 + (-1.6)  alpha_L * x - ln x * 5 / 2 + C"
    using alpha_L  _ C  _ x by (intro diff_mono add_mono mult_right_mono) auto
  also have "chi_L k  1" if "k  {..<30}" for k :: nat
    using that unfolding lessThan_nat_numeral pred_numeral_simps arith_simps
    by (elim insertE) (auto simp: chi_L_def)
  hence "alpha_L * x - ln x * 5 / 2 + C  ψ x"
    using primes_psi_lower_bound[of x] x by (simp add: C_def period_def)
  finally show ?thesis
    by (simp add: mult_ac)
qed

context balanced_chebyshev_multiset
begin

lemma psi_upper_bound:
  fixes x c C :: real
  defines "x0  Max ({3, 55 * c}  {3 * ¦l¦ |l. l ∈# L})"
  assumes x: "x  x0"
  assumes chi_nonneg: "n. n  {0..<period}  χL (real n)  0"
  assumes chi_ge1: "n. real n  {1..<c}  χL (real n)  1"
  assumes c: "c > 1" "c  period"
  assumes "αL  0"
  shows "ψ x  c / (c - 1) * αL * x + (3 * size L) / (4 * ln c) * ln x ^ 2 + ψ x0"
proof -
  have L_nonzero': "l  0" if "l ∈# L" for l
    using that L_nonzero by auto

  have chi_nonneg: "χL x  0" for x
    proof -
    have "χL x = χL (floor x mod period)"
      by (simp add: chi_L_mod chi_L_floor)
    also have "floor x mod period = real (nat (floor x mod period))"
      using period_pos by auto
    also have "χL   0"
      by (rule chi_nonneg) (use period_pos in auto simp: nat_less_iff)
    finally show ?thesis .
  qed

  have chi_ge1: "χL x  1" if "x  1" "x < c" for x
    proof -
    have "χL x = χL (floor x mod period)"
      by (simp add: chi_L_mod chi_L_floor)
    also have "floor x mod period = real (nat (floor x mod period))"
      using period_pos by auto
    also have "χL   1"
    proof (rule chi_ge1)
      have "real_of_int x < c"
        using that by linarith
      hence "real_of_int (x mod int period) < c"
        using that period_pos c by simp
      moreover have "1  x mod int period"
        by (use period_pos c that in auto simp: floor_less_iff)
      ultimately show "real (nat (x mod int period))  {1..<c}"
        by auto
    qed
    finally show ?thesis .
  qed

  have "finite {3 * l |l. l ∈# L}"
    by auto
  have x1: "x0  3" "x0  55 * c"
    unfolding x0_def by (rule Max_ge; simp)+
  have x2: "3 * ¦l¦  x0" if "l ∈# L" for l
    unfolding x0_def by (rule Max_ge) (use that in auto)

  define C where "C = 1/2 * (l∈#L. ln  ¦l¦) - size L"
  have *: "x  x * y" if "y  1" "x  0" for x y :: real
    using mult_left_mono[of 1 y x] that by simp

  have rec: "ψ x  ψ (x / c) + αL * x + ln x * size L / 2 - C" if x: "x  x0" for x :: real
  proof -
    have "x / c  x"
      using c using divide_left_mono[of 1 c x] x0  3 x by auto
    have "ψ x = ψ (x / c) + (d | d > 0  real d  {x/c<..x}. mangoldt d)"
      unfolding ψ_def sum_upto_def
      by (rule sum.union_disjoint' [symmetric])
         (use c x / c  x in auto)
    also have "(d | d > 0  real d  {x/c<..x}. mangoldt d) 
               (d | d > 0  real d  {x/c<..x}. mangoldt d * χL (x / d))"
      using c by (intro sum_mono * mangoldt_nonneg) (auto intro!: chi_ge1 simp: field_simps)
    also have "  (d | d > 0  real d  x. mangoldt d * χL (x / d))"
      by (intro sum_mono2) (auto intro!: mult_nonneg_nonneg mangoldt_nonneg chi_nonneg)
    also have " = ψL x"
      by (simp add: psi_L_def sum_upto_def)
    finally have "ψ x  ψ (x / c) + ψL x"
      by - simp_all

    have L: "3 * ¦real_of_int l¦  x" if "l ∈# L" for l
      using x2[OF that] x by linarith
  
    have "ψ x  ψ (x / c) + ψL x"
      by fact
    also have "ψL x = αL * x + psi_L_remainder x"
      using x0  3 x by (subst psi_L_eq) auto
    also have "¦psi_L_remainder x¦  ln x * size L / 2 - C"
      using psi_L_remainder_bounds[of x] x0  3 x L by (simp add: C_def)
    hence "psi_L_remainder x  ln x * size L / 2 - C"
      by linarith
    finally show "ψ x  ψ (x / c) + αL * x + ln x * size L / 2 - C"
      by (simp add: algebra_simps)
  qed

  define m where "m = nat log c (x / x0)"
  have "x > 0"
    using x x1 by simp

  have "ψ x  ψ x0 + (k<m. αL * x / c ^ k + ln (x / c ^ k) * size L / 2 - C)"
  proof -
    have "ψ x  ψ (x / c ^ m) + (k<m. αL * (x / c ^ k) + ln (x / c ^ k) * size L / 2 - C)"
    proof (rule primes_psi_lower_rec)
      fix x :: real assume "x  x0"
      thus "ψ x  ψ (x / c) + (αL * x + ln x * size L / 2 - C)"
        using rec[of x] by (simp add: algebra_simps)
    next
      have "c ^ m = c powr real m"
        using c by (simp add: powr_realpow)
      also have "  c powr (log c (x / x0) + 1)"
        using c x x0  3 by (intro powr_mono) (auto simp: m_def)
      also have " = c * x / x0"
        using c x x0  3 by (auto simp: powr_add)
      finally show "x0 * c ^ m  x * c"
        using x0  3 by (simp add: field_simps)
    qed (use x1 c in auto)
    also have "ψ (x / c ^ m)  ψ x0"
    proof (rule ψ_mono)
      have "x / x0 = c powr log c (x / x0)"
        using c x x0  3 by simp
      also have "  c powr m"
        unfolding m_def using c x0  3 x by (intro powr_mono) auto
      also have " = c ^ m"
        using c by (simp add: powr_realpow)
      finally show "x / c ^ m  x0"
        using x0  3 c by (simp add: field_simps)
    qed
    finally show ?thesis
      by simp
  qed
  also have " = ψ x0 + (k<m. αL * x / c ^ k + (ln x - k * ln c) * size L / 2 - C)"
    using x(1) x0  3 c by (simp add: ln_div ln_realpow)
  also have " = ψ x0 + αL * x * (k<m. 1 / c ^ k) + ln x * m * size L / 2 - real (k<m. k) * ln c * size L / 2 - C * m"
    by (simp add: sum_diff_distrib sum_subtractf sum.distrib sum_distrib_left sum_distrib_right algebra_simps diff_divide_distrib sum_divide_distrib)
  also have "(k<m. 1 / c ^ k) = (1 - (1 / c) ^ m) / (1 - 1 / c)"
    using sum_gp_strict[of "1/c" m] c by (simp add: field_simps)
  also have "  1 / (1 - 1 / c)"
    using c by (intro divide_right_mono) auto
  also have "1 / (1 - 1/c) = c / (c - 1)"
    using c by (simp add: field_simps)
  also have "(k<m. k) = real m * (real m - 1) / 2"
    by (induction m) (auto simp: field_simps)
  finally have "ψ x  ψ x0 + c / (c - 1) * αL * x +
                  ln x * m * size L / 2 -
                  real m * (real m - 1) / 2 * ln c * size L / 2 - C * m"
    using αL  0 x > 0 x1 by (simp add: mult_left_mono mult_right_mono mult_ac)
  also have " = ψ x0 + c / (c - 1) * αL * x + m/2 * (size L * (ln x - (real m - 1)/2 * ln c + 2) - (l∈#L. ln ¦l¦))"
    by (simp add: algebra_simps C_def)
  also have "m/2 * (size L * (ln x - (real m - 1)/2 * ln c + 2) - (l∈#L. ln ¦l¦)) 
             m/2 * (size L * (3/2 * ln x) - 0)"
  proof (intro mult_left_mono diff_mono)
    have "real m  log c (x / x0)"
      using c x0  3 x unfolding m_def by auto
    hence "ln x - (real m - 1)/2 * ln c + 2 
           ln x - (log c (x / x0) - 1)/2 * ln c + 2"
      using c by (intro diff_mono add_mono mult_right_mono divide_right_mono) auto
    also have " = (ln x + ln x0 + (ln c + 4)) / 2"
      using c x x0  3 by (simp add: log_def ln_div field_simps)
    also have "ln x0  ln x"
      using x x1 by simp
    also have "ln c + 4  ln x"
    proof -
      have "exp (4 :: real)  55"
        by (approximation 10)
      hence "exp 4 * c  55 * c"
        using c by (intro mult_right_mono) auto
      also have "55 * c  x0"
        by fact
      also have "  x"
        by fact
      finally have "exp (ln c + 4)  exp (ln x)"
        unfolding exp_add using c x1 x by (simp add: mult_ac)
      thus ?thesis
        by (simp only: exp_le_cancel_iff)
    qed
    also have "(ln x + ln x + ln x) / 2 = 3 / 2 * ln x"
      by simp
    finally show "ln x - (real m - 1) / 2 * ln c + 2  3 / 2 * ln x"
      by - simp
  qed (auto intro!: sum_mset_nonneg simp: L_nonzero' Ints_nonzero_abs_ge1)
  also have "m / 2 * (size L * (3/2 * ln x) - 0) = 3 / 4 * m * size L * ln x"
    by simp
  also have "  3 / 4 * (ln x / ln c) * size L * ln x"
  proof (intro mult_left_mono mult_right_mono)
    have "real m  log c (x / x0) + 1"
      unfolding m_def using c x x0  3 by auto
    also have " / 2 = (ln x / ln c + (1 - log c x0)) / 2"
      using x0  3 x  x0 c 
      by (simp add: log_def ln_div field_simps)
    also have "1 - log c x0  0"
      using x1 c by simp
    finally show "real m  ln x / ln c" by - simp_all
  qed (use x x1 in auto)
  also have " = (3 * size L) / (4 * ln c) * ln x ^ 2"
    by (simp add: power2_eq_square)
  finally show "ψ x  c / (c - 1) * αL * x + (3 * size L) / (4 * ln c) * ln x ^ 2 + ψ x0"
    by (simp add: algebra_simps)
qed

end


subsection ‹Final results›

theorem psi_lower_ge_9:
  assumes x: "x  41"
  shows   "ψ x  0.9 * x"
proof (cases "x  900")
  case False
  have "x{real 41..real 900}. primes_psi x  0.9 * x"
    by (rule check_psi_lower_correct[where prec = 16]) eval
  from bspec[OF this, of x] show ?thesis
    using assms False by simp
next
  case x: True
  define f :: "real  real"
    where "f = (λx. 0.02128 * x - 2.5 * ln x - 1.6)"
  have "0  f 900"
    unfolding f_def by (approximation 10)
  also have "f 900  f x"
    using x
  proof (rule DERIV_nonneg_imp_nondecreasing, goal_cases)
    case (1 t)
    have "(f has_real_derivative (0.02128 - 2.5 / t)) (at t)"
      unfolding f_def using 1 by (auto intro!: derivative_eq_intros)
    moreover have "0.02128 - 2.5 / t  0"
      using 1 by (auto simp: field_simps)
    ultimately show ?case
      by blast
  qed
  finally have "0.9 * x  0.9 * x + f x"
    by linarith
  also have " = 0.92128 * x - 2.5 * ln x - 1.6"
    by (simp add: f_def)
  also have "  ψ x"
    by (rule psi_lower_bound_precise) (use x in auto)
  finally show ?thesis .
qed

theorem primes_theta_ge_82:
  assumes "x  97"
  shows   "θ x  0.82 * x"
proof (cases "x  46000")
  case False
  have "x{real 97..real 46000}. θ x  0.82 * x"
    by (rule check_theta_lower_correct[where prec = 20]) eval
  from bspec[OF this, of x] show ?thesis
    using False assms by simp
next
  case True
  with assms have x: "x  46000"
    by auto
  define f :: "real  real"
    where "f = (λx. 0.10128 * x - 2.5 * ln x - 2 * ln x * sqrt x - 1.6)"
  have "0  f 46000"
    unfolding f_def by (approximation 30)
  also have "f 46000  f x"
    using x
  proof (rule DERIV_nonneg_imp_nondecreasing, goal_cases)
    case (1 t)
    define D where "D = 0.10128 - 2.5 / t - 2 * sqrt t / t - ln t / sqrt t"
    have deriv: "(f has_real_derivative D) (at t)"
      unfolding f_def
      by (rule derivative_eq_intros refl | use 1 in force)+
         (simp add: field_simps D_def)
    have "0.10128 - D = 2.5 / t + 2 / sqrt t + ln t / sqrt t"
      using 1 by (simp add: D_def field_simps del: div_add div_diff  div_mult_self1 div_mult_self2  div_mult_self3  div_mult_self4)
    also have "  2.5 / 46000 + 2 / 214 + ln t / sqrt t"
      using 1 by (intro add_mono) (auto simp: real_le_rsqrt)

    also have "ln t / sqrt t  ln 46000 / sqrt 46000"
      using 1(1)
    proof (rule DERIV_nonpos_imp_nonincreasing, goal_cases)
      case (1 u)
      have "((λt. ln t / sqrt t) has_real_derivative ((2 - ln u) / (2 * u * sqrt u))) (at u)"
        by (rule derivative_eq_intros refl | use 1 in force)+
           (use 1 in simp add: field_simps)
      moreover {
        have "2  ln (10::real)"
          by (approximation 30)
        also have "  ln u"
          using 1 by simp
        finally have "ln u  2" .
      }      
      hence "((2 - ln u) / (2 * u * sqrt u))  0"
        using 1 by (intro divide_nonpos_nonneg) auto
      ultimately show ?case
        by blast
    qed
    also have "  0.0501"
      by (approximation 30)
    also have "2.5 / 46000 + 2 / 214 + 0.0501  (0.10128 :: real)"
      by simp
    finally have "D  0"
      by simp
    with deriv show ?case by blast
  qed

  finally have "0.82 * x  0.82 * x + f x"
    by linarith
  also have " = 0.92128 * x - 2.5 * ln x - 2 * ln x * sqrt x - 1.6"
    by (simp add: f_def)
  also have "  0.92128 * x - 2.5 * ln x - 1.6 + θ x - ψ x"
    using ψ_minus_θ_bound[of x] x by simp
  also have "0.92128 * x - 2.5 * ln x - 1.6  ψ x"
    by (rule psi_lower_bound_precise) (use x in auto)
  finally show ?thesis by simp
qed

corollary primorial_ge_exp_82:
  assumes "x  97"
  shows   "primorial x  exp (0.82 * x)"
proof -
  have "primorial x = exp (θ x)"
    using ln_primorial[of x] primorial_pos[of x]
    by (metis exp_ln of_nat_0_less_iff)
  also have "  exp (0.82 * x)"
    using primes_theta_ge_82[OF assms] by simp
  finally show ?thesis .
qed


theorem primes_psi_le_111:
  assumes "x  0"
  shows   "ψ x  1.11 * x"
proof -
  have "x{real 0..real 146000}. primes_psi x  1.04 * x"
  proof (rule check_psi_upper_correct[where prec = 16])
    show "check_psi_upper (104 / 102) 16 0 146000"
      by eval
  qed auto
  hence initial: "primes_psi x  1.04 * x" if "x  {0..146000}" for x
    using that by auto

  show ?thesis
  proof (cases "x  146000")
    case False
    thus ?thesis
      using initial[of x] assms by simp
  next
    case x: True
    define L :: "int multiset" where "L = {#1, -2, -3, -5, 30#}"
    have [simp]: "set_mset L = {1, -2, -3, -5, 30}" "size L = 5"
      by (simp_all add: L_def)
    interpret balanced_chebyshev_multiset L
      by unfold_locales (auto simp: L_def)
    define x0 :: real where "x0 = Max ({3, 55 * 6}  {3 * ¦real_of_int l¦ |l. l ∈# L})"
  
    have x0: "x0 = 330"
    proof -
      have "x0 = Max ({3, 55 * 6}  {3 * ¦real_of_int l¦ |l. l ∈# L})"
        unfolding x0_def ..
      also have "{3 * ¦real_of_int l¦ |l. l ∈# L} = (λl. 3 * ¦of_int l¦) ` set_mset L"
        by blast
      finally show ?thesis
        by simp
    qed

    define f :: "real  real"
      where "f = (λt. 2.093 * ln t ^ 2 + 343.2 - 0.0044 * t)"
  
    have alpha_L: "alpha_L = ln 2 / 2 - (ln 30 / 30 - ln 5 / 5 - ln 3 / 3)"
      unfolding alpha_L_def by (simp add: L_def)
    have "alpha_L  0"
      unfolding alpha_L by (approximation 10)
    have period: "period = 30"
      by (simp add: period_def)
  
    have "ψ x  6 / (6 - 1) * alpha_L * x + (3 * size L) / (4 * ln 6) * ln x ^ 2 + ψ x0"
      unfolding x0_def
    proof (rule psi_upper_bound; (unfold period)?)
      show "chi_L (real n)  0" if "n  {0..<30}" for n
        unfolding chi_L_def
        using that unfolding atLeastLessThan_nat_numeral pred_numeral_simps arith_simps
        by (auto simp: L_def)
    next
      show "chi_L (real n)  1" if "real n  {1..<6}" for n
      proof -
        have "n  {1..<6}"
          using that by auto
        also have "{1..<6} = {1,2,3,4,5::nat}"
          by auto
        finally show ?thesis
          unfolding chi_L_def by (elim insertE) (auto simp: L_def)
      qed
    qed (use alpha_L  0 x in auto)
  
    also have " = 6/5 * alpha_L * x + 15 / (4 * ln 6) * (ln x)2 + ψ x0"
      by simp
    also have "ψ x0  343.2"
      using initial[of x0] by (simp add: x0)
    also have "6/5 * alpha_L  1.1056"
      unfolding alpha_L by (approximation 30)
    also have "15 / (4 * ln 6 :: real)  2.093"
      by (approximation 20)
    finally have "ψ x  1.1056 * x + 2.093 * ln x ^ 2 + 343.2"
      using x by - simp_all
    also have " = 1.11 * x + f x"
      by (simp add: algebra_simps f_def)
    also have "f x  f 146000"
      using x
    proof (rule DERIV_nonpos_imp_nonincreasing, goal_cases)
      case (1 t)
      define f' :: "real  real"
        where "f' = (λt. 4.186 * ln t / t - 0.0044)"
      have "(f has_field_derivative f' t) (at t)"
        using 1 unfolding f_def f'_def
        by (auto intro!: derivative_eq_intros)
      moreover {
        have "f' t  f' 146000"
          using 1(1)
        proof (rule DERIV_nonpos_imp_nonincreasing, goal_cases)
          case (1 u)
          have "(f' has_field_derivative (4.186 * (1 - ln u) / u ^ 2)) (at u)"
            using 1 unfolding f'_def
            by (auto intro!: derivative_eq_intros simp: field_simps power2_eq_square)
          moreover have "4.186 * (1 - ln u) / u ^ 2  0"
            using 1 exp_le by (auto intro!: divide_nonpos_nonneg simp: ln_ge_iff)
          ultimately show ?case
            by blast
        qed
        also have "  0"
          unfolding f'_def by (approximation 10)
        finally have "f' t  0" .
      }
      ultimately show ?case
        by blast
    qed
    also have "f 146000  0"
      unfolding f_def by (approximation 10)
    finally show ?thesis
      by - simp_all
  qed
qed

corollary primes_theta_le_111:
  assumes "x  0"
  shows   "θ x  1.11 * x"
  using primes_psi_le_111[OF assms] θ_le_ψ[of x]
  by linarith

text ‹
  As an easy corollary, we obtain Bertrand's postulate:
  For any real number $x > 1$, the interval $(x, 2x)$ contains
  at least one prime.
›
corollary bertrands_postulate:
  assumes "x > 1"
  shows   "p. prime p  real p  {x<..<2*x}"
proof (cases "x  7")
  case False
  consider "x  {1<..<2}" | "x  {2..<3}" | "x  {3..<5}" | "x  {5..<7}"
    using False assms by force
  thus ?thesis
  proof cases
    case 1
    thus ?thesis by (intro exI[of _ 2]; simp)
  next
    case 2
    thus ?thesis by (intro exI[of _ 3]; simp)
  next
    case 3
    thus ?thesis by (intro exI[of _ 5]; simp)
  next
    case 4
    thus ?thesis by (intro exI[of _ 7]; simp)
  qed
next
  case x: True
  have fin: "finite {p. prime p  real p  1.999 * x}"
    by (rule finite_subset[of _ "{..nat 2*x}"])
       (use x in auto simp: le_nat_iff le_floor_iff)

  have "θ (1.999 * x) > 1.11 * x"
  proof (cases "x  49")
    case False
    have "x{real 11..real 100}. θ x  0.556 * x"
      by (rule check_theta_lower_correct[where prec = 10]) eval
    from bspec[OF this, of "1.999*x"] show ?thesis
      using False x by simp
  next
    case True
    thus ?thesis
      using primes_theta_ge_82[of "1.999*x"] True by auto
  qed

  have "θ x  1.11 * x"
    by (rule primes_theta_le_111) (use x in auto)
  also have " < θ (1.999 * x)"
    by fact
  finally have "θ (1.999 * x) > θ x" .

  have "{p. prime p  real p  {x<..1.999*x}}  {}"
  proof
    assume eq: "{p. prime p  real p  {x<..1.999*x}} = {}"
    have "θ (1.999 * x) = θ x + (p | prime p  real p  {x<..1.999*x}. ln p)"
      unfolding primes_theta_def prime_sum_upto_def
      by (rule sum.union_disjoint' [symmetric]) (use fin in auto)
    also note eq
    finally show False
      using θ (1.999 * x) > θ x by simp
  qed
  thus ?thesis
    by auto
qed

unbundle no prime_counting_syntax

end