Theory Selberg_Asymptotic_Formula

(*
  File:    Selberg_Asymptotic_Formula.thy
  Author:  Manuel Eberl, TU München

  Selberg's asymptotic formula, which is an important ingredient in the elementary
  proof of the Prime Number Theorem by Selberg and Erdős.
*)
section ‹Selberg's asymptotic formula›
theory Selberg_Asymptotic_Formula
imports
  More_Dirichlet_Misc
  Prime_Number_Theorem.Prime_Counting_Functions
  Shapiro_Tauberian
  Euler_MacLaurin.Euler_MacLaurin_Landau
  Partial_Zeta_Bounds
begin

text ‹
  Following Apostol, we first show an inversion formula: Consider a function $f(x)$ for 
  $x\in\mathbb{R}_{{}>\,0}$. Define $g(x) := \ln x \cdot \sum_{n\leq x} f(x / n)$. Then:
  \[f(x) \ln x + \sum_{n\leq x} \Lambda(n) f(x/n) = \sum_{n\leq x} \mu(n) g(x/n)\]
›
(* 4.17 *)
locale selberg_inversion =
  fixes F G :: "real  'a :: {real_algebra_1, comm_ring_1}"
  defines "G  (λx. of_real (ln x) * sum_upto (λn. F (x / n)) x)"
begin  

lemma eq:
  assumes "x  1"
  shows "F x * of_real (ln x) + dirichlet_prod' mangoldt F x = dirichlet_prod' moebius_mu G x"
proof -
  have "F x * of_real (ln x) =
          dirichlet_prod' (λn. if n = 1 then 1 else 0) (λx. F x * of_real (ln x)) x"
    by (subst dirichlet_prod'_one_left) (use x  1 in auto)
  also have " = dirichlet_prod' (λn. d | d dvd n. moebius_mu d) (λx. F x * of_real (ln x)) x"
    by (intro dirichlet_prod'_cong refl, subst sum_moebius_mu_divisors') auto
  finally have eq1: "F x * of_real (ln x) = " .

  have eq2: "dirichlet_prod' mangoldt F x =
               dirichlet_prod' (dirichlet_prod moebius_mu (λn. of_real (ln (real n)))) F x"
  proof (intro dirichlet_prod'_cong refl)
    fix n :: nat assume n: "n > 0"
    thus "mangoldt n = dirichlet_prod moebius_mu (λn. of_real (ln (real n)) :: 'a) n"
      by (intro moebius_inversion mangoldt_sum [symmetric]) auto
  qed

  have "F x * of_real (ln x) + dirichlet_prod' mangoldt F x =
          sum_upto (λn. F (x / n) * (d | d dvd n.
            moebius_mu d * of_real (ln (x / n) + ln (n div d)))) x"
    unfolding eq1 eq2 unfolding dirichlet_prod'_def sum_upto_def
    by (simp add: algebra_simps sum.distrib dirichlet_prod_def sum_distrib_left sum_distrib_right)
  also have " = sum_upto (λn. F (x / n) * (d | d dvd n. moebius_mu d * of_real (ln (x / d)))) x"
    using x  1 by (intro sum_upto_cong refl arg_cong2[where f = "λx y. x * y"] sum.cong)
                     (auto elim!: dvdE simp: ln_div ln_mult)
  also have " = sum_upto (λn. d | d dvd n. moebius_mu d * of_real (ln (x / d)) * F (x / n)) x"
    by (simp add: sum_distrib_left sum_distrib_right mult_ac)
  also have " = ((n,d)(SIGMA n:{n. n > 0  real n  x}. {d. d dvd n}).
                    moebius_mu d * of_real (ln (x / d)) * F (x / n))"
    unfolding sum_upto_def by (subst sum.Sigma) (auto simp: case_prod_unfold)
  also have " = ((d,q)(SIGMA d:{d. d > 0  real d  x}. {q. q > 0  real q  x / d}).
                    moebius_mu d * of_real (ln (x / d)) * F (x / (q * d)))"
    by (rule sum.reindex_bij_witness[of _ "λ(d,q). (d * q, d)" "λ(n,d). (d, n div d)"])
       (auto simp: Real.real_of_nat_div field_simps dest: dvd_imp_le)
  also have " = sum_upto (λd. moebius_mu d * of_real (ln (x / d)) *
                                  sum_upto (λq. F (x / (q * d))) (x / d)) x"
    by (subst sum.Sigma [symmetric]) (auto simp: sum_upto_def sum_distrib_left)
  also have " = dirichlet_prod' moebius_mu G x"
    by (simp add: dirichlet_prod'_def G_def mult_ac)
  finally show ?thesis .
qed

end

text ‹
  We can now show Selberg's formula
  \[\psi(x)\ln x + \sum_{n\leq x} \Lambda(n)\psi(x/n) = 2x\ln x + O(x)\ .\]
›
(* 4.18 *)
theorem selberg_asymptotic_formula:
  includes prime_counting_syntax
  shows    "(λx. ψ x * ln x + dirichlet_prod' mangoldt ψ x) =o
               (λx. 2 * x * ln x) +o O(λx. x)"
proof -
  define C :: real where [simp]: "C = euler_mascheroni"
  define F2 :: "real  real" where [simp]: "F2 = (λx. x - C - 1)"
  define G1 where "G1 = (λx. ln x * sum_upto (λn. ψ (x / n)) x)"
  define G2 where "G2 = (λx. ln x * sum_upto (λn. F2 (x / n)) x)"

  interpret F1: selberg_inversion ψ G1
    by unfold_locales (simp_all add: G1_def)
  interpret F2: selberg_inversion F2 G2
    by unfold_locales (simp_all add: G2_def)

  have G1_bigo: "(λx. G1 x - (x * ln x ^ 2 - x * ln x))  O(λx. ln x ^ 2)"
  proof -
    have "(λx. ln x * (sum_upto (λn. ψ (x / n)) x - x * ln x + x))  O(λx. ln x * ln x)" 
      by (intro landau_o.big.mult_left sum_upto_ψ_x_over_n_asymptotics)
    thus ?thesis by (simp add: power2_eq_square G1_def algebra_simps)
  qed

  have G2_bigo: "(λx. G2 x - (x * ln x ^ 2 - x * ln x))  O(ln)"
  proof -
  define R1 :: "real  real" where "R1 = (λx. x * ln x * (harm (nat x) - (ln x + C)))"
  define R2 :: "real  real" where "R2 = (λx. (C + 1) * ln x * frac x)"
    have "(λx. G2 x - (x * ln x ^ 2 - x * ln x))  Θ(λx. R1 x + R2 x)"
    proof (intro bigthetaI_cong eventually_mono[OF eventually_ge_at_top[of 1]])
      fix x :: real assume x: "x  1"
      have "G2 x = x * ln x * sum_upto (λn. 1 / n) x - (C + 1) * x * ln x"
        using x by (simp add: G2_def sum_upto_altdef sum_subtractf
                              sum_distrib_left sum_distrib_right algebra_simps)
      also have "sum_upto (λn. 1 / n) x = harm (nat x)"
        using x unfolding sum_upto_def harm_def
        by (intro sum.cong) (auto simp: field_simps le_nat_iff le_floor_iff)
      also have "x * ln x * harm (nat x) - (C + 1) * x * ln x =
                   x * ln x ^ 2 - x * ln x + R1 x + R2 x"
        by (simp add: R1_def R2_def algebra_simps frac_def power2_eq_square)
      finally show "G2 x - (x * ln x ^ 2 - x * ln x) = R1 x + R2 x" by simp
    qed
    also have "(λx. R1 x + R2 x)  O(ln)"
    proof (intro sum_in_bigo)
      have "(λx::real. ln x - ln (nat x))  O(λx. ln x - ln (x - 1))"
      proof (intro bigoI[of _ 1] eventually_mono[OF eventually_ge_at_top[of 2]])
        fix x :: real assume x: "x  2"
        thus "norm (ln x - ln (nat x))  1 * norm (ln x - ln (x - 1))" by auto
      qed
      also have "(λx::real. ln x - ln (x - 1))  O(λx. 1 / x)" by real_asymp
      finally have bigo_ln_floor: "(λx::real. ln x - ln (nat x))  O(λx. 1 / x)" .
  
      have "(λx. harm (nat x) - (ln (nat x) + C))  O(λx. 1 / nat x)"
        unfolding C_def using harm_expansion_bigo_simple2
        by (rule landau_o.big.compose) (* TODO: real_asymp should be able to do this *)
           (auto intro!: filterlim_compose[OF filterlim_nat_sequentially filterlim_floor_sequentially])
      also have "(λx. 1 / nat x)  O(λx. 1 / x)" by real_asymp
      finally have "(λx. harm (nat x) - (ln (nat x) + C) - (ln x - ln (nat x)))
                        O(λx. 1 / x)" by (rule sum_in_bigo[OF _bigo_ln_floor])
      hence "(λx. harm (nat x) - (ln x + C))  O(λx. 1 / x)" by (simp add: algebra_simps)
      hence "(λx. x * ln x * (harm (nat x) - (ln x + C)))  O(λx. x * ln x * (1 / x))"
        by (intro landau_o.big.mult_left)
      thus "R1  O(ln)" by (simp add: landau_divide_simps R1_def)
    next
      have "R2  O(λx. 1 * ln x * 1)" (* TODO: real_asymp? *)
        unfolding R2_def by (intro landau_o.big.mult landau_o.big_refl) real_asymp+
      thus "R2  O(ln)" by (simp add: R2_def)
    qed
    finally show "(λx. G2 x - (x * (ln x)2 - x * ln x))  O(ln)" .
  qed
  hence G2_bigo': "(λx. G2 x - (x * (ln x)2 - x * ln x))  O(λx. ln x ^ 2)"
    by (rule landau_o.big.trans) real_asymp+

  ― ‹Now things become a bit hairy. In order to show that the `Big-O' bound is actually
      valid for all x ≥ 1›, we need to show that G1 x - G2 x› is bounded on any compact
      interval starting at 1.›
  have "c>0. x1. ¦G1 x - G2 x¦  c * ¦sqrt x¦"
  proof (rule bigoE_bounded_real_fun)
    have "(λx. G1 x - G2 x)  O(λx. ln x ^ 2)"
      using sum_in_bigo(2)[OF G1_bigo G2_bigo'] by simp
    also have "(λx::real. ln x ^ 2)  O(sqrt)" by real_asymp
    finally show "(λx. G1 x - G2 x)  O(sqrt)" .
  next
    fix x :: real assume "x  1"
    thus "¦sqrt x¦  1" by simp
  next
    fix b :: real assume b: "b  1"
    show "bounded ((λx. G1 x - G2 x) ` {1..b})"
    proof (rule boundedI, safe)
      fix x assume x: "x  {1..b}"
      have "¦G1 x - G2 x¦ = ¦ln x * sum_upto (λn. ψ (x / n) - F2 (x / n)) x¦"
        by (simp add: G1_def G2_def sum_upto_def sum_distrib_left ring_distribs sum_subtractf)
      also have " = ln x * ¦sum_upto (λn. ψ (x / n) - F2 (x / n)) x¦"
        using x b by (simp add: abs_mult)
      also have "¦sum_upto (λn. ψ (x / n) - F2 (x / n)) x¦ 
                   sum_upto (λn. ¦ψ (x / n) - F2 (x / n)¦) x"
        unfolding sum_upto_def by (rule sum_abs)
      also have "  sum_upto (λn. ψ x + (x + C + 1)) x"
        unfolding sum_upto_def
      proof (intro sum_mono)
        fix n assume n: "n  {n. n > 0  real n  x}"
        hence le: "x / n  x / 1" by (intro divide_left_mono) auto
        thus "¦ψ (x / n) - F2 (x / n)¦  ψ x + (x + C + 1)"
          unfolding F2_def using euler_mascheroni_pos x le ψ_nonneg ψ_mono[of "x / n" x]
          by (intro order.trans[OF abs_triangle_ineq]
                    order.trans[OF abs_triangle_ineq4] add_mono) auto
      qed
      also have " = (ψ x + (x + C + 1)) * x"
        using x by (simp add: sum_upto_altdef)
      also have "ln x * ((ψ x + (x + C + 1)) * real_of_int x) 
                   ln b * ((ψ b + (b + C + 1)) * real_of_int b)"
        using euler_mascheroni_pos x
        by (intro mult_mono add_mono order.refl ψ_mono add_nonneg_nonneg mult_nonneg_nonneg
                     ψ_nonneg) (auto intro: floor_mono)
      finally show "norm (G1 x - G2 x)  ln b * ((ψ b + (b + C + 1)) * real_of_int b)"
        using x by (simp add: mult_left_mono)
    qed
  qed auto
  then obtain A where A: "A > 0" "x. x  1  ¦G1 x - G2 x¦  A * sqrt x" by auto

  ― ‹The rest of the proof now consists merely of combining some asymptotic estimates.›
  have "(λx. (ψ x - F2 x) * ln x + sum_upto (λn. mangoldt n * (ψ (x / n) - F2 (x / n))) x)
           Θ(λx. sum_upto (λn. moebius_mu n * (G1 (x / n) - G2 (x / n))) x)"
  proof (intro bigthetaI_cong eventually_mono[OF eventually_ge_at_top[of 1]])
    fix x :: real assume x: "x  1"
    have "(ψ x - F2 x) * ln x + sum_upto (λn. mangoldt n * (ψ (x / n) - F2 (x / n))) x =
            (ψ x * of_real (ln x) + dirichlet_prod' mangoldt ψ x) -
            (F2 x * of_real (ln x) + dirichlet_prod' mangoldt F2 x)"
      by (simp add: algebra_simps dirichlet_prod'_def sum_upto_def sum_subtractf sum.distrib)
    also have " = sum_upto (λn. moebius_mu n * (G1 (x / n) - G2 (x / n))) x"
      unfolding F1.eq[OF x] F2.eq[OF x]
      by (simp add: dirichlet_prod'_def sum_upto_def sum_subtractf sum.distrib algebra_simps)
    finally show "(ψ x - F2 x) * ln x + sum_upto (λn. mangoldt n * (ψ (x/n) - F2 (x/n))) x = " .
  qed
  also have "(λx. sum_upto (λn. moebius_mu n * (G1 (x / n) - G2 (x / n))) x) 
               O(λx. A * sqrt x * sum_upto (λx. x powr (-1/2)) x)"
  proof (intro bigoI eventually_mono[OF eventually_ge_at_top[of 1]])
    fix x :: real assume x: "x  1"
    have "¦sum_upto (λn. moebius_mu n * (G1 (x / n) - G2 (x / n))) x¦ 
            sum_upto (λn. ¦moebius_mu n * (G1 (x / n) - G2 (x / n))¦) x"
      unfolding sum_upto_def by (rule sum_abs)
    also have "  sum_upto (λn. 1 * (A * sqrt (x / n))) x"
      unfolding sum_upto_def abs_mult by (intro A sum_mono mult_mono) (auto simp: moebius_mu_def)
    also have " = A * sqrt x * sum_upto (λx. x powr (-1/2)) x"
      using x by (simp add: sum_upto_def powr_minus powr_half_sqrt sum_distrib_left
                            sum_distrib_right real_sqrt_divide field_simps)
    also have "  ¦A * sqrt x * sum_upto (λx. x powr (-1/2)) x¦" by simp
    finally show "norm (sum_upto (λn. moebius_mu n * (G1 (x / n) - G2 (x / n))) x) 
                    1 * norm (A * sqrt x * sum_upto (λx. x powr (-1/2)) x)" by simp
  qed
  also have "(λx. A * sqrt x * sum_upto (λx. x powr (-1/2)) x)  O(λx. 1 * sqrt x * x powr (1/2))"
    using zeta_partial_sum_le_pos_bigo[of "1 / 2"]
    by (intro landau_o.big.mult ) (auto simp: max_def)
  also have "(λx::real. 1 * sqrt x * x powr (1/2))  O(λx. x)"
    by real_asymp
  finally have bigo: "(λx. (ψ x - F2 x) * ln x + sum_upto (λn. mangoldt n * (ψ (x/n) - F2 (x/n))) x)
                         O(λx. x)" (is "?h  _") .

  let ?R = "λx. sum_upto (λn. mangoldt n / n) x"
  let ?lhs = "λx. ψ x * ln x + dirichlet_prod' mangoldt ψ x"

  note bigo
  also have "?h = (λx. ?lhs x - (x * ln x - (C + 1) * (ln x + ψ x)) - x * ?R x)"
    by (rule ext) (simp add: algebra_simps dirichlet_prod'_def  sum_distrib_right ψ_def
                             sum_upto_def sum_subtractf sum.distrib sum_distrib_left)
  finally have "(λx. ?lhs x - (x * ln x - (C + 1) * (ln x + ψ x)) - x * ?R x + x * (?R x - ln x))
                     O(λx. x)" (is "?h'  _")
  proof (rule sum_in_bigo)
    have "(λx. x * (sum_upto (λn. mangoldt n / real n) x - ln x))  O(λx. x * 1)"
      by (intro landau_o.big.mult_left ψ.asymptotics)
    thus "(λx. x * (sum_upto (λn. mangoldt n / real n) x - ln x))  O(λx. x)" by simp
  qed
  also have "?h' = (λx. ?lhs x - (2 * x * ln x - (C + 1) * (ln x + ψ x)))"
    by (simp add: fun_eq_iff algebra_simps)
  finally have "(λx. ?lhs x - (2*x*ln x - (C+1) * (ln x + ψ x)) - (C+1) * (ln x + ψ x))  O(λx. x)"
  proof (rule sum_in_bigo)
    have "(λx. ln x + ψ x)  O(λx. x)"
      by (intro sum_in_bigo bigthetaD1[OF ψ.bigtheta]) real_asymp+
    thus "(λx. (C + 1) * (ln x + ψ x))  O(λx. x)" by simp
  qed
  also have "(λx. ?lhs x - (2*x*ln x - (C+1) * (ln x + ψ x)) - (C+1) * (ln x + ψ x)) =
               (λx. ?lhs x - 2 * x * ln x)" by (simp add: algebra_simps)
  finally show ?thesis
    by (subst set_minus_plus [symmetric]) (simp_all add: fun_diff_def algebra_simps)
qed

end