Theory Summatory_Divisor_Sigma_Bounds

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

  Asymptotic expansions for the sums over the divisor function σ_x.
*)
section ‹The asymptotics of the summatory divisor $\sigma$ function›
theory Summatory_Divisor_Sigma_Bounds
  imports Partial_Zeta_Bounds More_Dirichlet_Misc
begin

text ‹
  In this section, we analyse the asymptotic behaviour of the summatory divisor functions
  $\sum_{n\leq x} \sigma_\alpha(n)$ for real α›. This essentially tells us what the average value
  of these functions is for large x›.

  The case α = 0› is not treated here since $\sigma_0$ is simply the divisor function,
  for which precise asymptotics are already available in the AFP.
›

subsection ‹Case 1: $\alpha = 1$›

text ‹
  If α = 1›, $\sigma_\alpha(n)$ is simply the sum of all divisors of n›. Here, the asymptotics is
  \[\sum_{n\leq x} \sigma_1(n) = \frac{\pi^2}{12} x^2 + O(x\ln x)\ .\]
›
theorem summatory_divisor_sum_asymptotics:
  "sum_upto divisor_sum =o (λx. pi2 / 12 * x ^ 2) +o O(λx. x * ln x)"
proof -
  define ζ where "ζ = Re (zeta 2)"
  define R1 where "R1 = (λx. sum_upto real x - x2 / 2)"
  define R2 where "R2 = (λx. sum_upto (λd. 1 / d2) x - (ζ - 1 / x))"
  obtain c1 where c1: "c1 > 0" "x. x  1  ¦R1 x¦  c1 * x"
    using zeta_partial_sum_le_neg[of 1] by (auto simp: R1_def)
  obtain c2 where c2: "c2 > 0" "x. x  1  ¦R2 x¦  c2 / x2"
    using zeta_partial_sum_le_pos[of 2]
     by (auto simp: ζ_def R2_def powr_minus field_simps
              simp del: div_mult_self3 div_mult_self4 div_mult_self2 div_mult_self1)

  have le: "¦sum_upto divisor_sum x - ζ / 2 * x2¦  c2 / 2 + x / 2 + c1 * x * (ln x + 1)"
    if x: "x  1" for x
  proof -
    have div_le: "real (a div b)  x" if "a  x" for a b :: nat
      by (rule order.trans[OF _ that(1)]) auto
 
    have "real (sum_upto divisor_sum x) = sum_upto (dirichlet_prod real (λ_. 1)) x"
      by (simp add: divisor_sigma_conv_dirichlet_prod [abs_def] 
                    sum_upto_def divisor_sigma_1_left [symmetric])
    also have " = sum_upto (λn. d | d dvd n. real d) x"
      by (simp add: dirichlet_prod_def)
    also have " = ((n, d)  (SIGMA n:{n. n > 0  real n  x}. {d. d dvd n}). real d)"
      unfolding sum_upto_def by (subst sum.Sigma) auto
    also have " = ((d, q)  (SIGMA d:{d. d > 0  real d  x}. {q. q > 0  real q  x / d}). real q)"
      by (rule sum.reindex_bij_witness[of _ "λ(d, q). (d * q, q)" "λ(n, d). (n div d, d)"])
         (use div_le in auto simp: field_simps)
    also have " = sum_upto (λd. sum_upto real (x / d)) x"
      by (simp add: sum_upto_def sum.Sigma)
    also have " = x2 * sum_upto (λd. 1 / d2) x / 2 + sum_upto (λd. R1 (x / d)) x"
      by (simp add: R1_def sum_upto_def sum.distrib sum_subtractf sum_divide_distrib
                    power_divide sum_distrib_left)
    also have "sum_upto (λd. 1 / d2) x = ζ - 1 / x + R2 x"
      by (simp add: R2_def)
    finally have eq: "real (sum_upto divisor_sum x) =
                        x2 * (ζ - 1 / x + R2 x) / 2 + sum_upto (λd. R1 (x / real d)) x" .
  
    have "real (sum_upto divisor_sum x) - ζ / 2 * x2 =
            x2 / 2 * R2 x - x / 2 + sum_upto (λd. R1 (x / real d)) x" using x
      by (subst eq)
         (simp add: field_simps power2_eq_square del: div_diff div_add
               del: div_mult_self3 div_mult_self4 div_mult_self2 div_mult_self1)
    also have "¦¦  c2 / 2 + x / 2 + c1 * x * (ln x + 1)"
    proof (rule order.trans[OF abs_triangle_ineq] order.trans[OF abs_triangle_ineq4] add_mono)+
      have "¦x2 / 2 * R2 x¦ = x2 / 2 * ¦R2 x¦"
        using x by (simp add: abs_mult)
      also have "  x2 / 2 * (c2 / x2)"
        using x by (intro mult_left_mono c2) auto
      finally show "¦x2 / 2 * R2 x¦  c2 / 2"
        using x by simp
    next
      have "¦sum_upto (λd. R1 (x / real d)) x¦  sum_upto (λd. ¦R1 (x / real d)¦) x"
        unfolding sum_upto_def by (rule sum_abs)
      also have "  sum_upto (λd. c1 * (x / real d)) x"
        unfolding sum_upto_def by (intro sum_mono c1) auto
      also have " = c1 * x * sum_upto (λd. 1 / real d) x"
        by (simp add: sum_upto_def sum_distrib_left)
      also have "sum_upto (λd. 1 / real d) x = harm (nat x)"
        unfolding sum_upto_altdef harm_def by (intro sum.cong) (auto simp: field_simps)
      also have "  ln (nat x) + 1"
        by (rule harm_le) (use x in auto simp: le_nat_iff)
      also have "ln (nat x)  ln x" using x by simp
      finally show "¦sum_upto (λd. R1 (x / real d)) x¦  c1 * x * (ln x + 1)"
        using c1(1) x by simp
    qed (use x in auto)
    finally show "¦sum_upto divisor_sum x - ζ / 2 * x2¦  c2 / 2 + x / 2 + c1 * x * (ln x + 1)" .
  qed

  have "eventually (λx. ¦sum_upto divisor_sum x - ζ / 2 * x2¦ 
                            c2 / 2 + x / 2 + c1 * x * (ln x + 1)) at_top"
    using eventually_ge_at_top[of 1] by eventually_elim (use le in auto)
  hence "eventually (λx. ¦sum_upto divisor_sum x - ζ / 2 * x2¦ 
                              ¦c2 / 2 + x / 2 + c1 * x * (ln x + 1)¦) at_top"
    by eventually_elim linarith
  hence "(λx. sum_upto divisor_sum x - ζ / 2 * x2)  O(λx. c2 / 2 + x / 2 + c1 * x * (ln x + 1))"
    by (intro landau_o.bigI[of 1]) auto
  also have "(λx. c2 / 2 + x / 2 + c1 * x * (ln x + 1))  O(λx. x * ln x)"
    by real_asymp
  finally show ?thesis
    by (subst set_minus_plus [symmetric])
       (simp_all add: fun_diff_def algebra_simps ζ_def zeta_even_numeral)
qed


subsection ‹Case 2: $\alpha > 0$, $\alpha \neq 1$›

(* TODO: Unify *)
text ‹
  Next, we consider the case $\alpha > 0$ and $\alpha\neq 1$. We then have:
  \[\sum_{n\leq x} \sigma_\alpha(n) = \frac{\zeta(\alpha + 1)}{\alpha + 1} x^{\alpha + 1} +
       O\left(x^{\text{max}(1,\alpha)}\right)\]
›
theorem summatory_divisor_sigma_asymptotics_pos:
  fixes α :: real
  assumes α: "α > 0" "α  1"
  defines "ζ  Re (zeta (α + 1))"
  shows  "sum_upto (divisor_sigma α) =o
            (λx. ζ / (α + 1) * x powr (α + 1)) +o O(λx. x powr max 1 α)"
proof -
  define R1 where "R1 = (λx. sum_upto (λd. real d powr α) x - x powr (α + 1) / (α + 1))"
  define R2 where "R2 = (λx. sum_upto (λd. d powr (-α - 1)) x - (ζ - x powr -α / α))"
  define R3 where "R3 = (λx. sum_upto (λd. d powr -α) x - x powr (1 - α) / (1 - α))"
  obtain c1 where c1: "c1 > 0" "x. x  1  ¦R1 x¦  c1 * x powr α"
    using zeta_partial_sum_le_neg[of α] α by (auto simp: R1_def add_ac)
  obtain c2 where c2: "c2 > 0" "x. x  1  ¦R2 x¦  c2 * x powr (-α-1)"
    using zeta_partial_sum_le_pos[of "α + 1"] α by (auto simp: ζ_def R2_def)
  obtain c3 where c3: "c3 > 0" "x. x  1  ¦R3 x¦  c3"
    using zeta_partial_sum_le_pos'[of "α"] α by (auto simp: R3_def)
  define ub :: "real  real" where
    "ub = (λx. x / (α * (α + 1)) + c2 / (α + 1) + c1 * (1 / (1 - α) * x + c3 * x powr α))"

  have le: "¦sum_upto (divisor_sigma α) x - ζ / (α + 1) * x powr (α + 1)¦  ub x"
    if x: "x  1" for x
  proof -
    have div_le: "real (a div b)  x" if "a  x" for a b :: nat
      by (rule order.trans[OF _ that(1)]) auto
 
    have "sum_upto (divisor_sigma α) x =
            sum_upto (dirichlet_prod (λn. real n powr α) (λ_. 1)) x"
      by (simp add: divisor_sigma_conv_dirichlet_prod [abs_def] 
                    sum_upto_def divisor_sigma_1_left [symmetric])
    also have " = sum_upto (λn. d | d dvd n. real d powr α) x"
      by (simp add: dirichlet_prod_def)
    also have " = ((n, d)  (SIGMA n:{n. n > 0  real n  x}. {d. d dvd n}). real d powr α)"
      unfolding sum_upto_def by (subst sum.Sigma) auto
    also have " = ((d, q)  (SIGMA d:{d. d > 0  real d  x}. {q. q > 0  real q  x / d}). real q powr α)"
      by (rule sum.reindex_bij_witness[of _ "λ(d, q). (d * q, q)" "λ(n, d). (n div d, d)"])
         (use div_le in auto simp: field_simps)
    also have " = sum_upto (λd. sum_upto (λq. q powr α) (x / d)) x"
      by (simp add: sum_upto_def sum.Sigma)
    also have " = x powr (α + 1) * sum_upto (λd. 1 / d powr (α + 1)) x / (α + 1) +
                    sum_upto (λd. R1 (x / d)) x"
      by (simp add: R1_def sum_upto_def sum.distrib sum_subtractf sum_divide_distrib
                    powr_divide sum_distrib_left)
    also have "sum_upto (λd. 1 / d powr (α + 1)) x = ζ - x powr -α / α + R2 x"
      by (simp add: R2_def powr_minus field_simps powr_diff powr_add)
    finally have eq: "sum_upto (divisor_sigma α) x =
       x powr (α + 1) * (ζ - x powr -α / α + R2 x) / (α + 1) + sum_upto (λd. R1 (x / d)) x" .
  
    have "sum_upto (divisor_sigma α) x - ζ / (α + 1) * x powr (α + 1) =
            -x / (α * (α + 1)) + x powr (α + 1) / (α + 1) * R2 x + sum_upto (λd. R1 (x / d)) x"
      using x α
      by (subst eq, simp add: divide_simps
                         del: div_mult_self3 div_mult_self4 div_mult_self2 div_mult_self1)
         (simp add: field_simps power2_eq_square powr_add powr_minus del: div_diff div_add
               del: div_mult_self3 div_mult_self4 div_mult_self2 div_mult_self1)
    also have "¦¦  ub x" unfolding ub_def
    proof (rule order.trans[OF abs_triangle_ineq] order.trans[OF abs_triangle_ineq4] add_mono)+
      have "¦x powr (α + 1) / (α + 1) * R2 x¦ = x powr (α + 1) / (α + 1) * ¦R2 x¦"
        using x α by (simp add: abs_mult)
      also have "  x powr (α + 1) / (α + 1) * (c2 * x powr (-α-1))"
        using x α by (intro mult_left_mono c2) auto
      also have " = c2 / (α + 1)"
        using α x by (simp add: field_simps powr_diff powr_minus powr_add)
      finally show "¦x powr (α + 1) / (α + 1) * R2 x¦  c2 / (α + 1)" .
    next
      have "¦sum_upto (λd. R1 (x / real d)) x¦  sum_upto (λd. ¦R1 (x / real d)¦) x"
        unfolding sum_upto_def by (rule sum_abs)
      also have "  sum_upto (λd. c1 * (x / real d) powr α) x"
        unfolding sum_upto_def by (intro sum_mono c1) auto
      also have " = c1 * x powr α * sum_upto (λd. 1 / real d powr α) x"
        by (simp add: sum_upto_def sum_distrib_left powr_divide)
      also have "sum_upto (λd. 1 / real d powr α) x = x powr (1-α) / (1-α) + R3 x"
        using x by (simp add: R3_def powr_minus field_simps)
      also have "c1 * x powr α * (x powr (1 - α) / (1 - α) + R3 x) =
                   c1 / (1 - α) * x + c1 * x powr α * R3 x"
        using x by (simp add: powr_diff divide_simps
                    del: div_mult_self3 div_mult_self4 div_mult_self2 div_mult_self1)
                   (simp add: field_simps)
      also have "c1 * x powr α * R3 x  c1 * x powr α * c3"
        using x c1(1) c3(2)[of x] by (intro mult_left_mono) auto
      finally show "¦sum_upto (λd. R1 (x / d)) x¦  c1 * (1 / (1 - α) * x + c3 * x powr α)"
        by (simp add: field_simps)
    qed (use α x in simp_all)
    finally show "¦sum_upto (divisor_sigma α) x - ζ / (α + 1) * x powr (α + 1)¦  ub x" .
  qed

  have "eventually (λx. ¦sum_upto (divisor_sigma α) x - ζ / (α+1) * x powr (α+1)¦  ub x) at_top"
    using eventually_ge_at_top[of 1] by eventually_elim (use le in auto)
  hence "eventually (λx. ¦sum_upto (divisor_sigma α) x - ζ/(α+1) * x powr (α+1)¦  ¦ub x¦) at_top"
    by eventually_elim linarith
  hence "(λx. sum_upto (divisor_sigma α) x - ζ/(α+1) * x powr (α+1))  O(ub)"
    by (intro landau_o.bigI[of 1]) auto
  also have "ub  O(λx. x powr max 1 α)"
    using α unfolding ub_def by (cases "α  1"; real_asymp)
  finally show ?thesis
    by (subst set_minus_plus [symmetric])
       (simp_all add: fun_diff_def algebra_simps ζ_def zeta_even_numeral)
qed


subsection ‹Case 3: $\alpha < 0$›

text ‹
  Last, we consider the case of a negative exponent. We have for $\alpha > 0$:
  \[\sum_{n\leq x} \sigma_{-\alpha}(n) = \zeta(\alpha + 1)x + O(R(x))\]
  where $R(x) = \ln x$ if $\alpha = 1$ and $R(x) = x^{\text{max}(0, 1-\alpha)}$ otherwise.
›
theorem summatory_divisor_sigma_asymptotics_neg:
  fixes α :: real
  assumes α: "α > 0"
  defines "δ  max 0 (1 - α)"
  defines "ζ  Re (zeta (α + 1))"
  shows  "sum_upto (divisor_sigma (-α)) =o (if α = 1 then (λx. pi2/6 * x) +o O(ln)
                                                     else (λx. ζ * x) +o O(λx. x powr δ))"
proof -
  define Ra where "Ra = (λx. -sum_upto (λd. d powr (-α) * frac (x / d)) x)"
  define R1 where "R1 = (λx. sum_upto (λd. real d powr (-α)) x - (x powr (1 - α) / (1 - α) + ζ))"
  define R2 where "R2 = (λx. sum_upto (λd. d powr (-α - 1)) x - (ζ - x powr -α / α))"
  define R3 where "R3 = (λx. sum_upto (λd. d powr -α) x - x powr (1 - α) / (1 - α))"
  obtain c2 where c2: "c2 > 0" "x. x  1  ¦R2 x¦  c2 * x powr (-α-1)"
    using zeta_partial_sum_le_pos[of "α + 1"] α by (auto simp: ζ_def R2_def)
  define ub :: "real  real" where
    "ub = (λx. x powr (1 - α) / α + c2 * x powr - α + ¦Ra x¦)"

  have le: "¦sum_upto (divisor_sigma (-α)) x - ζ * x¦  ub x"
    if x: "x  1" for x
  proof -
    have div_le: "real (a div b)  x" if "a  x" for a b :: nat
      by (rule order.trans[OF _ that(1)]) auto
 
    have "sum_upto (divisor_sigma (-α)) x =
            sum_upto (dirichlet_prod (λn. real n powr (-α)) (λ_. 1)) x"
      by (simp add: divisor_sigma_conv_dirichlet_prod [abs_def] 
                    sum_upto_def divisor_sigma_1_left [symmetric])
    also have " = sum_upto (λn. d | d dvd n. real d powr (-α)) x"
      by (simp add: dirichlet_prod_def)
    also have " = ((n, d)  (SIGMA n:{n. n > 0  real n  x}. {d. d dvd n}). real d powr (-α))"
      unfolding sum_upto_def by (subst sum.Sigma) auto
    also have " = ((d, q)  (SIGMA d:{d. d > 0  real d  x}. {q. q > 0  real q  x / d}).
                         real d powr (-α))"
      by (rule sum.reindex_bij_witness[of _ "λ(d, q). (d * q, d)" "λ(n, d). (d, n div d)"])
         (use div_le in auto simp: field_simps dest: dvd_imp_le)
    also have " = sum_upto (λd. sum_upto (λq. d powr (-α)) (x / d)) x"
      by (simp add: sum_upto_def sum.Sigma [symmetric])
    also have " = sum_upto (λd. d powr (-α) * x / d) x"
      using x by (simp add: sum_upto_altdef mult_ac)
    also have " = x * sum_upto (λd. d powr (-α) / d) x + Ra x"
      by (simp add: frac_def sum_distrib_left sum_distrib_right
                    sum_subtractf sum_upto_def algebra_simps Ra_def)
    also have "sum_upto (λd. d powr (-α) / d) x = sum_upto (λd. d powr (-α - 1)) x"
      by (simp add: powr_diff powr_minus powr_add field_simps)
    also have " = ζ - x powr - α / α + R2 x"
      by (simp add: R2_def)
    finally have "sum_upto (divisor_sigma (-α)) x - ζ * x = -(x powr (1 - α) / α) + x * R2 x + Ra x"
      using x α by (simp add: powr_diff powr_minus field_simps)

    also have "¦¦  x powr (1 - α) / α + c2 * x powr -α + ¦Ra x¦"
    proof (rule order.trans[OF abs_triangle_ineq] order.trans[OF abs_triangle_ineq4] add_mono)+
      from x have "¦x * R2 x¦  x * ¦R2 x¦"
        by (simp add: abs_mult)
      also from x have "  x * (c2 * x powr (-α - 1))"
        by (intro mult_left_mono c2) auto
      also have " = c2 * x powr -α"
        using x by (simp add: field_simps powr_minus powr_diff)
      finally show "¦x * R2 x¦  " .
    qed (use x α in auto)
    finally show "¦sum_upto (divisor_sigma (- α)) x - ζ * x¦  ub x"
      by (simp add: ub_def)
  qed

  have "eventually (λx. ¦sum_upto (divisor_sigma (-α)) x - ζ * x¦  ub x) at_top"
    using eventually_ge_at_top[of 1] by eventually_elim (use le in auto)
  hence "eventually (λx. ¦sum_upto (divisor_sigma (-α)) x - ζ * x¦  ¦ub x¦) at_top"
    by eventually_elim linarith
  hence bigo: "(λx. sum_upto (divisor_sigma (-α)) x - ζ * x)  O(ub)"
    by (intro landau_o.bigI[of 1]) auto

  define ub' :: "real  real" where "ub' = sum_upto (λn. real n powr - α)"
  have "¦Ra x¦  ¦ub' x¦" if "x  1" for x
  proof -
    have "¦Ra x¦  sum_upto (λn. ¦real n powr -α * frac (x / n)¦) x"
      unfolding Ra_def abs_minus sum_upto_def by (rule sum_abs)
    also have "  sum_upto (λn. real n powr -α * 1) x"
      unfolding abs_mult sum_upto_def
      by (intro sum_mono mult_mono) (auto intro: less_imp_le[OF frac_lt_1])
    finally show ?thesis by (simp add: ub'_def)
  qed
  hence "Ra  O(ub')"
    by (intro bigoI[of _ 1] eventually_mono[OF eventually_ge_at_top[of 1]]) auto
  also have "ub'  O(λx. if α = 1 then ln x else x powr δ)"
  proof (cases "α = 1")
    case [simp]: True
    have "sum_upto (λn. 1 / n)  O(ln)"
      by (intro asymp_equiv_imp_bigo harm_asymp_equiv)
    thus ?thesis by (simp add: ub'_def powr_minus field_simps)
  next
    case False
    have "sum_upto (λn. real n powr - α)  O(λx. x powr δ)"
      using assms False unfolding δ_def by (intro zeta_partial_sum_pos_bigtheta bigthetaD1)
    thus ?thesis
      using zeta_partial_sum_neg_asymp_equiv[of α] α False by (simp add: ub'_def)
  qed
  finally have Ra_bigo: "Ra  " .

  show ?thesis
  proof (cases "α = 1")
    case [simp]: True
    with Ra_bigo have Ra: "(λx. ¦Ra x¦)  O(ln)" by simp
    note bigo
    also have "ub  O(λx. ln x)"
      unfolding ub_def by (intro sum_in_bigo Ra) real_asymp+
    finally have "sum_upto (divisor_sigma (-α)) =o (λx. (pi2 / 6) * x) +o O(ln)"
      by (subst set_minus_plus [symmetric])
         (simp_all add: fun_diff_def algebra_simps ζ_def zeta_even_numeral)
    thus ?thesis by (simp only: True refl if_True)
  next
    case False
    with Ra_bigo have Ra: "(λx. ¦Ra x¦)  O(λx. x powr δ)" by simp
    have *: "(λx. x powr (1 - α) / α)  O(λx. x powr δ)"
            "(λx. c2 * x powr - α)  O(λx. x powr δ)"
      unfolding δ_def using α False by (cases "α > 1"; real_asymp)+

    note bigo
    also have "ub  O(λx. x powr δ)"
      unfolding ub_def using α False by (intro sum_in_bigo Ra *)
    finally have "sum_upto (divisor_sigma (-α)) =o (λx. ζ * x) +o O(λx. x powr δ)"
      by (subst set_minus_plus [symmetric])
         (simp_all add: fun_diff_def algebra_simps ζ_def zeta_even_numeral)
    thus ?thesis by (simp only: False refl if_False)
  qed
qed

end