Theory Arithmetic_Summatory_Asymptotics

(*
    File:      Arithmetic_Summatory_Asymptotics.thy
    Author:    Manuel Eberl, TU München
*)
section ‹Asymptotics of summatory arithmetic functions›
theory Arithmetic_Summatory_Asymptotics
  imports 
    Euler_MacLaurin.Euler_MacLaurin_Landau
    Arithmetic_Summatory 
    Dirichlet_Series_Analysis
    Landau_Symbols.Landau_More
begin

subsection ‹Auxiliary bounds›

lemma sum_inverse_squares_tail_bound:
  assumes "d > 0"
  shows   "summable (λn. 1 / (real (Suc n) + d) ^ 2)"
          "(n. 1 / (real (Suc n) + d) ^ 2)  1 / d"
proof -
  show *: "summable (λn. 1 / (real (Suc n) + d) ^ 2)"
  proof (rule summable_comparison_test, intro allI exI impI)
    fix n :: nat
    from assms show "norm (1 / (real (Suc n) + d) ^ 2)  1 / real (Suc n) ^ 2"
      unfolding norm_divide norm_one norm_power
      by (intro divide_left_mono power_mono) simp_all
  qed (insert inverse_squares_sums, simp add: sums_iff)
  show "(n. 1 / (real (Suc n) + d) ^ 2)  1 / d"
  proof (rule sums_le)
    fix n have "1 / (real (Suc n) + d) ^ 2  1 / ((real n + d) * (real (Suc n) + d))"
      unfolding power2_eq_square using assms
      by (intro divide_left_mono mult_mono mult_pos_pos add_nonneg_pos) simp_all
    also have " = 1 / (real n + d) - 1 / (real (Suc n) + d)"
      using assms by (simp add: divide_simps)
    finally show "1 / (real (Suc n) + d)2  1 / (real n + d) - 1 / (real (Suc n) + d)" .
  next
    show "(λn. 1 / (real (Suc n) + d)2) sums (n. 1 / (real (Suc n) + d)2)"
      using * by (simp add: sums_iff)
  next
    have "(λn. 1 / (real n + d) - 1 / (real (Suc n) + d)) sums (1 / (real 0 + d) - 0)"
      by (intro telescope_sums' real_tendsto_divide_at_top[OF tendsto_const],
          subst add.commute, rule filterlim_tendsto_add_at_top[OF tendsto_const 
            filterlim_real_sequentially])
    thus "(λn. 1 / (real n + d) - 1 / (real (Suc n) + d)) sums (1 / d)" by simp
  qed
qed

lemma moebius_sum_tail_bound:
  assumes "d > 0"
  shows   "abs (n. moebius_mu (Suc n + d) / real (Suc n + d)^2)  1 / d" (is "abs ?S  _")
proof -
  have *: "summable (λn. 1 / (real (Suc n + d))2)"
    by (insert sum_inverse_squares_tail_bound(1)[of "real d"] assms, simp_all add: add_ac)
  have **: "summable (λn. abs (moebius_mu (Suc n + d) / real (Suc n + d)^2))"
  proof (rule summable_comparison_test, intro exI allI impI)
    fix n :: nat
    show "norm (¦moebius_mu (Suc n + d) / (real (Suc n + d))^2¦) 
            1 / (real (Suc n + d))^2" 
      unfolding real_norm_def abs_abs abs_divide power_abs abs_of_nat
      by (intro divide_right_mono abs_moebius_mu_le) simp_all
  qed (insert *)   
  from ** have "abs ?S  (n. abs (moebius_mu (Suc n + d) / real (Suc n + d)^2))"
    by (rule summable_rabs)
  also have "  (n. 1 / (real (Suc n) + d) ^ 2)"
  proof (intro suminf_le allI)
    fix n :: nat
    show "abs (moebius_mu (Suc n + d) / (real (Suc n + d))^2)  1 / (real (Suc n) + real d)^2"
      unfolding abs_divide abs_of_nat power_abs of_nat_add [symmetric]
      by (intro divide_right_mono abs_moebius_mu_le) simp_all
  qed (insert * **, simp_all add: add_ac)
  also from assms have "  1 / d" by (intro sum_inverse_squares_tail_bound) simp_all
  finally show ?thesis .
qed

lemma sum_upto_inverse_bound:
  "sum_upto (λi. 1 / real i) x  0"
  "eventually (λx. sum_upto (λi. 1 / real i) x  ln x + 13 / 22) at_top"
proof -
  show "sum_upto (λi. 1 / real i) x  0"
    by (simp add: sum_upto_def sum_nonneg)
  from order_tendstoD(2)[OF euler_mascheroni_LIMSEQ euler_mascheroni_less_13_over_22]
  obtain N where N: "n. n  N  harm n - ln (real n) < 13 / 22"
    unfolding eventually_at_top_linorder by blast
  show "eventually (λx. sum_upto (λi. 1 / real i) x  ln x + 13 / 22) at_top"
    using eventually_ge_at_top[of "max (real N) 1"]
  proof eventually_elim
    case (elim x)
    have "sum_upto (λi. 1 / real i) x = (i{0<..nat x}. 1 / real i)"
      by (simp add: sum_upto_altdef)
    also have " = harm (nat x)"
      unfolding harm_def by (intro sum.cong refl) (auto simp: field_simps)
    also have "  ln (real (nat x)) + 13 / 22"
      using N[of "nat x"] elim by (auto simp: le_nat_iff le_floor_iff)
    also have "ln (real (nat x))  ln x" using elim by (subst ln_le_cancel_iff) auto
    finally show ?case by - simp
  qed
qed

lemma sum_upto_inverse_bigo: "sum_upto (λi. 1 / real i)  O(λx. ln x)"
proof -
  have "eventually (λx. norm (sum_upto (λi. 1 / real i) x)  1 * norm (ln x + 13/22)) at_top"
    using eventually_ge_at_top[of "1::real"] sum_upto_inverse_bound(2)
    by eventually_elim (insert sum_upto_inverse_bound(1), simp_all)
  hence "sum_upto (λi. 1 / real i)  O(λx. ln x + 13/22)"
    by (rule bigoI)
  also have "(λx::real. ln x + 13/22)  O(λx. ln x)" by simp
  finally show ?thesis .
qed

lemma
  defines "G  (λx::real. (n. moebius_mu (n + Suc (nat x)) / (n + Suc (nat x))^2) :: real)"
  shows   moebius_sum_tail_bound': "t. t  2  ¦G t¦  1 / (t - 1)"
    and   moebius_sum_tail_bigo:   "G  O(λt. 1 / t)"
proof -
  show "¦G t¦  1 / (t - 1)" if t: "t  2" for t
  proof -
    from t have "¦G t¦  1 / real (nat t)"
      unfolding G_def using moebius_sum_tail_bound[of "nat t"] by simp
    also have "t  1 + real_of_int t" by linarith
    hence "1 / real (nat t)  1 / (t - 1)" using t by (simp add: field_simps)
    finally show ?thesis .
  qed
  hence "G  O(λt. 1 / (t - 1))" 
    by (intro bigoI[of _ 1] eventually_mono[OF eventually_ge_at_top[of "2::real"]]) auto
  also have "(λt::real. 1 / (t - 1))  Θ(λt. 1 / t)" by simp
  finally show "G  O(λt. 1 / t)" .
qed


subsection ‹Summatory totient function›

theorem summatory_totient_asymptotics:
  "(λx. sum_upto (λn. real (totient n)) x - 3 / pi2 * x2)  O(λx. x * ln x)"
proof -
  define H where "H = (λx. of_int (floor x) * (of_int (floor x) + 1) / 2 - x ^ 2 / 2 :: real)"
  define H' where "H' = (λx. sum_upto (λn. moebius_mu n * H (x / real n)) x)"
  have H: "sum_upto real x = x^2/2 + H x" if "x  0" for x
    using that by (simp add: sum_upto_real H_def)
  define G where "G = (λx::real. (n. moebius_mu (n + Suc (nat x)) / (n + Suc (nat x))^2))"
  
  have H_bound: "¦H t¦  t / 2" if "t  0" for t
  proof -
    have "H t - t / 2 = (-(t - of_int (floor t))) * (floor t + t + 1) / 2"
      by (simp add: H_def field_simps power2_eq_square)
    also have "  0" using that by (intro mult_nonpos_nonneg divide_nonpos_nonneg) simp_all
    finally have "H t  t / 2" by simp
    have "-H t - t / 2 = (t - of_int (floor t) - 1) * (of_int (floor t) + t) / 2"
      by (simp add: H_def field_simps power2_eq_square)
    also have "  0" using that
      by (intro divide_nonpos_nonneg mult_nonpos_nonneg) ((simp; fail) | linarith)+
    finally have "-H t  t / 2" by simp
    with H t  t / 2 show "¦H t¦  t / 2" by simp
  qed
  
  have H'_bound: "¦H' t¦  t / 2 * sum_upto (λi. 1 / real i) t" if "t  0" for t
  proof -
    have "¦H' t¦  (i | 0 < i  real i  t. ¦moebius_mu i * H (t / real i)¦)"
      unfolding H'_def sum_upto_def by (rule sum_abs)
    also have "  (i | 0 < i  real i  t. 1 * ((t / real i) / 2))"
      unfolding abs_mult using that
      by (intro sum_mono mult_mono abs_moebius_mu_le H_bound) simp_all
    also have " = t / 2 * sum_upto (λi. 1 / real i) t"
      by (simp add: sum_upto_def sum_distrib_left sum_distrib_right mult_ac)
    finally show ?thesis .
  qed
  hence "H'  O(λt. t * sum_upto (λi. 1 / real i) t)"
    using sum_upto_inverse_bound(1)
    by (intro bigoI[of _ "1/2"] eventually_mono[OF eventually_ge_at_top[of "0::real"]]) 
       (auto elim!: eventually_mono simp: abs_mult)
  also have "(λt. t * sum_upto (λi. 1 / real i) t)  O(λt. t * ln t)"
    by (intro landau_o.big.mult sum_upto_inverse_bigo) simp_all
  finally have H'_bigo: "H'  O(λx. x * ln x)" .
  
  {
    fix x :: real assume x: "x  0"
    have "sum_upto (λn. real (totient n)) x = sum_upto (λn. of_int (int (totient n))) x"
      by simp
    also have " = sum_upto (λn. moebius_mu n * sum_upto real (x / real n)) x"
      by (subst totient_conv_moebius_mu) (simp add: sum_upto_dirichlet_prod of_int_dirichlet_prod)
    also have " = sum_upto (λn. moebius_mu n * ((x / real n) ^ 2 / 2 + H (x / real n))) x" using x
      by (intro sum_upto_cong) (simp_all add: H)
    also have " = x^2 / 2 * sum_upto (λn. moebius_mu n / real n ^ 2) x + H' x"
      by (simp add: sum_upto_def H'_def sum.distrib ring_distribs
                    sum_distrib_left sum_distrib_right power_divide mult_ac)
    also have "sum_upto (λn. moebius_mu n / real n ^ 2) x = 
                 (n{..<Suc (nat x)}. moebius_mu n / real n ^ 2)" 
      unfolding sum_upto_altdef by (intro sum.mono_neutral_cong_left refl) auto
    also have " = 6 / pi ^ 2 - G x"
      using sums_split_initial_segment[OF moebius_over_square_sums, of "Suc (nat x)"]
      by (auto simp: sums_iff algebra_simps G_def)
    finally have "sum_upto (λn. real (totient n)) x = 3 / pi2 * x2 - x2 / 2 * G x + H' x"
      by (simp add: algebra_simps)
  }
  hence "(λx. sum_upto (λn. real (totient n)) x - 3 / pi^2 * x^2)  
            Θ(λx. (-(x^2) / 2) * G x + H' x)"
    by (intro bigthetaI_cong eventually_mono[OF eventually_ge_at_top[of "0::real"]]) 
       (auto elim!: eventually_mono)
  also have "(λx. (-(x^2) / 2) * G x + H' x)  O(λx. x * ln x)"
  proof (intro sum_in_bigo H'_bigo)
    have "(λx. (- (x^2) / 2) * G x)  O(λx. x^2 * (1 / x))"
      using moebius_sum_tail_bigo [folded G_def] by (intro landau_o.big.mult) simp_all
    also have "(λx::real. x^2 * (1 / x))  O(λx. x * ln x)" by simp
    finally show "(λx. (- (x^2) / 2) * G x)  O(λx. x * ln x)" .
  qed
  finally show ?thesis .
qed

theorem summatory_totient_asymptotics':
  "(λx. sum_upto (λn. real (totient n)) x) =o (λx. 3 / pi2 * x2) +o O(λx. x * ln x)"
  using summatory_totient_asymptotics
  by (subst set_minus_plus [symmetric]) (simp_all add: fun_diff_def)

theorem summatory_totient_asymptotics'':
  "sum_upto (λn. real (totient n)) ∼[at_top] (λx. 3 / pi2 * x2)"
proof -
  have "(λx. sum_upto (λn. real (totient n)) x - 3 / pi2 * x2)  O(λx. x * ln x)"
    by (rule summatory_totient_asymptotics)
  also have "(λx. x * ln x)  o(λx. 3 / pi ^ 2 * x ^ 2)" by simp
  finally show ?thesis by (simp add: asymp_equiv_altdef)
qed


subsection ‹Asymptotic distribution of squarefree numbers›

lemma le_sqrt_iff: "x  0  x  sqrt y  x^2  y"
  using real_sqrt_le_iff[of "x^2" y] by (simp del: real_sqrt_le_iff)

theorem squarefree_asymptotics: "(λx. card {n. real n  x  squarefree n} - 6 / pi2 * x)  O(sqrt)"
proof -
  define f :: "nat  real" where "f = (λn. if n = 0 then 0 else 1)"
  define g :: "nat  real" where "g = dirichlet_prod (ind squarefree) moebius_mu"
  
  interpret g: multiplicative_function g unfolding g_def
    by (intro multiplicative_dirichlet_prod squarefree.multiplicative_function_axioms 
              moebius_mu.multiplicative_function_axioms)
  interpret g: multiplicative_function' g "λp k. if k = 2 then -1 else 0" "λ_. 0"
  proof
    interpret g': multiplicative_dirichlet_prod' "ind squarefree" moebius_mu
     "λp k. if 1 < k then 0 else 1" "λp k. if k = 1 then - 1 else 0" "λ_. 1" "λ_. - 1"
    by (intro multiplicative_dirichlet_prod'.intro squarefree.multiplicative_function'_axioms 
              moebius_mu.multiplicative_function'_axioms)
    fix p k :: nat assume "prime p" "k > 0"
    hence "g (p ^ k) = (i{0<..<k}. (if Suc 0 < i then 0 else 1) *
                           (if k - i = Suc 0 then - 1 else 0))"
      by (auto simp: g'.prime_power g_def)
    also have " = (i{0<..<k}. (if k = 2 then -1 else 0))" 
      by (intro sum.cong refl) auto
    also from k > 0 have " = (if k = 2 then -1 else 0)" by simp
    finally show "g (p ^ k) = " .
  qed simp_all
  have mult_g_square: "multiplicative_function (λn. g (n ^ 2))"
    by standard (simp_all add: power_mult_distrib g.mult_coprime)
    
  have g_square: "g (m ^ 2) = moebius_mu m" for m
    using mult_g_square moebius_mu.multiplicative_function_axioms
  proof (rule multiplicative_function_eqI)
    fix p k :: nat assume *: "prime p" "k > 0"
    have "g ((p ^ k) ^ 2) = g (p ^ (2 * k))" by (simp add: power_mult [symmetric] mult_ac)
    also from * have " = (if k = 1 then -1 else 0)" by (simp add: g.prime_power)
    also from * have " = moebius_mu (p ^ k)" by (simp add: moebius_mu.prime_power)
    finally show "g ((p ^ k) ^ 2) = moebius_mu (p ^ k)" .
  qed

  have g_nonsquare: "g m = 0" if "¬is_square m" for m
  proof (cases "m = 0")
    case False
    from that False obtain p where p: "prime p" "odd (multiplicity p m)"
      using is_nth_power_conv_multiplicity_nat[of 2 m] by auto
    from p have "multiplicity p m  2" by auto
    moreover from p have "p  prime_factors m" 
      by (auto simp: prime_factors_multiplicity intro!: Nat.gr0I)
    ultimately have "(pprime_factors m. if multiplicity p m = 2 then - 1 else 0 :: real) = 0"
      (is "?P = _") by auto
    also have "?P = g m" using False by (subst g.prod_prime_factors') auto
    finally show ?thesis .
  qed auto
    
  have abs_g_le: "abs (g m)  1" for m
    by (cases "is_square m")
       (auto simp: g_square g_nonsquare abs_moebius_mu_le elim!: is_nth_powerE)
    
  have fds_g: "fds g = fds_ind squarefree * fds moebius_mu"
    by (rule fds_eqI) (simp add: g_def fds_nth_mult)
  have "fds g * fds_zeta = fds_ind squarefree * (fds_zeta * fds moebius_mu)"
    by (simp add: fds_g mult_ac)
  also have "fds_zeta * fds moebius_mu = (1 :: real fds)"
    by (rule fds_zeta_times_moebius_mu)
  finally have *: "fds_ind squarefree = fds g * fds_zeta" by simp
  have ind_squarefree: "ind squarefree = dirichlet_prod g f"
  proof
    fix n :: nat
    from * show "ind squarefree n = dirichlet_prod g f n"
      by (cases "n = 0") (simp_all add: fds_eq_iff fds_nth_mult f_def)
  qed
    
  define H :: "real  real" 
    where "H = (λx. sum_upto (λm. g (m^2) * (real_of_int x / real (m2) - x / real (m^2))) (sqrt x))"
  define J where "J = (λx::real. (n. moebius_mu (n + Suc (nat x)) / (n + Suc (nat x))^2))"
      
  have "eventually (λx. norm (H x)  1 * norm (sqrt x)) at_top" 
    using eventually_ge_at_top[of "0::real"]
  proof eventually_elim
    case (elim x)
    have "abs (H x)  sum_upto (λm. abs (g (m^2) * (real_of_int x / real (m2) - 
                         x / real (m^2)))) (sqrt x)" (is "_  ?S") unfolding H_def sum_upto_def
      by (rule sum_abs)
    also have "x / (real m)2 - real_of_int x / (real m)2  1" for m by linarith
    hence "?S  sum_upto (λm. 1 * 1) (sqrt x)" unfolding abs_mult sum_upto_def
      by (intro sum_mono mult_mono abs_g_le) simp_all
    also have " = of_int sqrt x" using elim by (simp add: sum_upto_altdef)
    also have "  sqrt x" by linarith
    finally show ?case using elim by simp
  qed
  hence H_bigo: "H  O(λx. sqrt x)" by (rule bigoI)
  
  let ?A = "λx. card {n. real n  x  squarefree n}"
  have "eventually (λx. ?A x - 6 / pi2 * x = (-x) * J (sqrt x) + H x) at_top"
    using eventually_ge_at_top[of "0::real"]
  proof eventually_elim
    fix x :: real assume x: "x  0"
    have "{n. real n  x  squarefree n} = {n. n > 0  real n  x  squarefree n}" 
      by (auto intro!: Nat.gr0I)
    also have "card  = sum_upto (ind squarefree :: nat  real) x"
      by (rule sum_upto_ind [symmetric])
    also have " = sum_upto (λd. g d * sum_upto f (x / real d)) x" (is "_ = ?S")
      unfolding ind_squarefree by (rule sum_upto_dirichlet_prod)
    also have "sum f {0<..nat x / real i} = of_int x / real i" if "i > 0" for i
      using x by (simp add: f_def)
    hence "?S = sum_upto (λd. g d * of_int x / real d) x"
      unfolding sum_upto_altdef by (intro sum.cong refl) simp_all
    also have " = sum_upto (λm. g (m ^ 2) * of_int x / real (m ^ 2)) (sqrt x)"
      unfolding sum_upto_def 
    proof (intro sum.reindex_bij_betw_not_neutral [symmetric])
      show "bij_betw power2 ({i. 0 < i  real i  sqrt x} - {}) 
              ({i. 0 < i  real i  x} - {i{0<..nat x}. ¬is_square i})"
        by (auto simp: bij_betw_def inj_on_def power_eq_iff_eq_base le_sqrt_iff 
                       is_nth_power_def le_nat_iff le_floor_iff)
    qed (auto simp: g_nonsquare)
    also have " = x * sum_upto (λm. g (m ^ 2) / real m ^ 2) (sqrt x) + H x"
      by (simp add: H_def sum_upto_def sum.distrib ring_distribs sum_subtractf 
                    sum_distrib_left sum_distrib_right mult_ac)
    also have "sum_upto (λm. g (m ^ 2) / real m ^ 2) (sqrt x) = 
                 sum_upto (λm. moebius_mu m / real m ^ 2) (sqrt x)"
      unfolding sum_upto_altdef by (intro sum.cong refl) (simp_all add: g_square)
    also have "sum_upto (λm. moebius_mu m / (real m)2) (sqrt x) =
                 (m<Suc (nat sqrt x). moebius_mu m / (real m) ^ 2)"
      unfolding sum_upto_altdef by (intro sum.mono_neutral_cong_left) auto
    also have " = (6 / pi^2 - J (sqrt x))"
      using sums_split_initial_segment[OF moebius_over_square_sums, of "Suc (nat sqrt x)"]
      by (auto simp: sums_iff algebra_simps J_def sum_upto_altdef)
    finally show "?A x - 6 / pi2 * x = (-x) * J (sqrt x) + H x"
      by (simp add: algebra_simps)
  qed
  hence "(λx. ?A x - 6 / pi2 * x)  Θ(λx. (-x) * J (sqrt x) + H x)"
    by (rule bigthetaI_cong)
  also have "(λx. (-x) * J (sqrt x) + H x)  O(λx. sqrt x)"
  proof (intro sum_in_bigo H_bigo)
    have "(λx. J (sqrt x))  O(λx. 1 / sqrt x)" unfolding J_def
        using moebius_sum_tail_bigo sqrt_at_top by (rule landau_o.big.compose)
    hence "(λx. (-x) * J (sqrt x))  O(λx. x * (1 / sqrt x))"
      by (intro landau_o.big.mult) simp_all
    also have "(λx::real. x * (1 / sqrt x))  Θ(λx. sqrt x)"
      by (intro bigthetaI_cong eventually_mono[OF eventually_gt_at_top[of "0::real"]]) 
         (auto simp: field_simps)
    finally show "(λx. (-x) * J (sqrt x))  O(λx. sqrt x)" .
  qed
  finally show ?thesis .
qed

theorem squarefree_asymptotics':
  "(λx. card {n. real n  x  squarefree n}) =o (λx. 6 / pi2 * x) +o O(λx. sqrt x)"
  using squarefree_asymptotics
  by (subst set_minus_plus [symmetric]) (simp_all add: fun_diff_def)

theorem squarefree_asymptotics'':
  "(λx. card {n. real n  x  squarefree n}) ∼[at_top] (λx. 6 / pi2 * x)"
proof -
  have "(λx. card {n. real n  x  squarefree n} - 6 / pi2 * x)  O(λx. sqrt x)"
    by (rule squarefree_asymptotics)
  also have "(sqrt :: real  real)  Θ(λx. x powr (1/2))" 
    by (intro bigthetaI_cong eventually_mono[OF eventually_ge_at_top[of "0::real"]]) 
       (auto simp: powr_half_sqrt)
  also have "(λx::real. x powr (1/2))  o(λx. 6 / pi ^ 2 * x)" by simp
  finally show ?thesis by (simp add: asymp_equiv_altdef)
qed


subsection ‹The hyperbola method›

lemma hyperbola_method_bigo:
  fixes f g :: "nat  'a :: real_normed_field"
  assumes "(λx. sum_upto (λn. f n * sum_upto g (x / real n)) (sqrt x) - R x)  O(b)"
  assumes "(λx. sum_upto (λn. sum_upto f (x / real n) * g n) (sqrt x) - S x)  O(b)"
  assumes "(λx. sum_upto f (sqrt x) * sum_upto g (sqrt x) - T x)  O(b)"
  shows   "(λx. sum_upto (dirichlet_prod f g) x - (R x + S x - T x))  O(b)"
proof -
  let ?A = "λx. (sum_upto (λn. f n * sum_upto g (x / real n)) (sqrt x) - R x) +
                (sum_upto (λn. sum_upto f (x / real n) * g n) (sqrt x) - S x) +
                (-(sum_upto f (sqrt x) * sum_upto g (sqrt x) - T x))"
  have "(λx. sum_upto (dirichlet_prod f g) x - (R x + S x - T x))  Θ(?A)"
    by (intro bigthetaI_cong eventually_mono[OF eventually_ge_at_top[of "0::real"]])
       (auto simp: hyperbola_method_sqrt)
  also from assms have "?A  O(b)"
    by (intro sum_in_bigo(1)) (simp_all only: landau_o.big.uminus_in_iff)
  finally show ?thesis .
qed

lemma frac_le_1: "frac x  1"
  unfolding frac_def by linarith

lemma ln_minus_ln_floor_bound:
  assumes "x  2"
  shows   "ln x - ln (floor x)  {0..<1 / (x - 1)}"
proof -
  from assms have "ln (floor x)  ln (x - 1)" by (subst ln_le_cancel_iff) simp_all
  hence "ln x - ln (floor x)  ln ((x - 1) + 1) - ln (x - 1)" by simp
  also from assms have " < 1 / (x - 1)" by (intro ln_diff_le_inverse) simp_all
  finally have "ln x - ln (floor x) < 1 / (x - 1)" by simp
  moreover from assms have "ln x  ln (of_int x)" by (subst ln_le_cancel_iff) simp_all
  ultimately show ?thesis by simp
qed
  
lemma ln_minus_ln_floor_bigo:
  "(λx::real. ln x - ln (floor x))  O(λx. 1 / x)"
proof -
  have "eventually (λx. norm (ln x - ln (floor x))  1 * norm (1 / (x - 1))) at_top"
    using eventually_ge_at_top[of "2::real"]
  proof eventually_elim
    case (elim x)
    with ln_minus_ln_floor_bound[OF this] show ?case by auto
  qed
  hence "(λx::real. ln x - ln (floor x))  O(λx. 1 / (x - 1))" by (rule bigoI) 
  also have "(λx::real. 1 / (x - 1))  O(λx. 1 / x)" by simp
  finally show ?thesis .
qed

lemma divisor_count_asymptotics_aux:
  "(λx. sum_upto (λn. sum_upto (λ_. 1) (x / real n)) (sqrt x) -
                    (x * ln x / 2 + euler_mascheroni * x))  O(sqrt)"
proof -
  define R where "R = (λx. i{0<..nat sqrt x}. frac (x / real i))"
  define S where "S = (λx. ln (real (nat sqrt x)) - ln x / 2)"
  have R_bound: "R x  {0..sqrt x}" if x: "x  0" for x
  proof -
    have "R x  (i{0<..nat sqrt x}. 1)" unfolding R_def by (intro sum_mono frac_le_1)
    also from x have " = of_int sqrt x" by simp
    also have "  sqrt x" by simp
    finally have "R x  sqrt x" .
    moreover have "R x  0" unfolding R_def by (intro sum_nonneg) simp_all
    ultimately show ?thesis by simp
  qed
  have R_bound': "norm (R x)  1 * norm (sqrt x)" if "x  0" for x 
    using R_bound[OF that] that by simp
  have R_bigo: "R  O(sqrt)" using eventually_ge_at_top[of "0::real"]
    by (intro bigoI[of _ 1], elim eventually_mono) (rule R_bound')
  
  have "eventually (λx. sum_upto (λn. sum_upto (λ_. 1 :: real) (x / real n)) (sqrt x) =
                          x * harm (nat sqrt x) - R x) at_top"
    using eventually_ge_at_top[of "0 :: real"]
  proof eventually_elim
    case (elim x)
    have "sum_upto (λn. sum_upto (λ_. 1 :: real) (x / real n)) (sqrt x) = 
            (i{0<..nat sqrt x}. of_int x / real i)" using elim
      by (simp add: sum_upto_altdef)
    also have " = x * (i{0<..nat sqrt x}. 1 / real i) - R x"
      by (simp add: sum_subtractf frac_def R_def sum_distrib_left)
    also have "{0<..nat sqrt x} = {1..nat sqrt x}" by auto
    also have "(i. 1 / real i) = harm (nat sqrt x)" by (simp add: harm_def divide_simps)
    finally show ?case .
  qed
  hence "(λx. sum_upto (λn. sum_upto (λ_. 1 :: real) (x / real n)) (sqrt x) -
                (x * ln x / 2 + euler_mascheroni * x))  
         Θ(λx. x * (harm (nat sqrt x) - (ln (nat sqrt x) + euler_mascheroni)) - R x + x * S x)" 
   (is "_  Θ(?A)")
   by (intro bigthetaI_cong) (elim eventually_mono, simp_all add: algebra_simps S_def)
  also have "?A  O(sqrt)"
  proof (intro sum_in_bigo)
    have "(λx. - S x)  Θ(λx. ln (sqrt x) - ln (of_int sqrt x))"
      by (intro bigthetaI_cong eventually_mono [OF eventually_ge_at_top[of "1::real"]]) 
         (auto simp: S_def ln_sqrt)
    also have "(λx. ln (sqrt x) - ln (of_int sqrt x))  O(λx. 1 / sqrt x)"
      by (rule landau_o.big.compose[OF ln_minus_ln_floor_bigo sqrt_at_top])
    finally have "(λx. x * S x)  O(λx. x * (1 / sqrt x))" by (intro landau_o.big.mult) simp_all
    also have "(λx::real. x * (1 / sqrt x))  Θ(λx. sqrt x)"
      by (intro bigthetaI_cong eventually_mono [OF eventually_gt_at_top[of "0::real"]]) 
         (auto simp: field_simps)
    finally show "(λx. x * S x)  O(sqrt)" .
  next
    let ?f = "λx::real. harm (nat sqrt x) - (ln (real (nat sqrt x)) + euler_mascheroni)"
    have "?f  O(λx. 1 / real (nat sqrt x))"
    proof (rule landau_o.big.compose[of _ _ _ "λx. nat sqrt x"])
      show "filterlim (λx::real. nat sqrt x) at_top at_top"
        by (intro filterlim_compose[OF filterlim_nat_sequentially]
                  filterlim_compose[OF filterlim_floor_sequentially] sqrt_at_top)
    next
      show "(λa. harm a - (ln (real a) + euler_mascheroni))  O(λa. 1 / real a)"
        by (rule harm_expansion_bigo_simple2)
    qed
    also have "(λx. 1 / real (nat sqrt x))  O(λx. 1 / (sqrt x - 1))"
    proof (rule bigoI[of _ 1], use eventually_ge_at_top[of 2] in eventually_elim)
      case (elim x)
      have "sqrt x  1 + real_of_int sqrt x" by linarith
      with elim show ?case by (simp add: field_simps)
    qed
    also have "(λx::real. 1 / (sqrt x - 1))  O(λx. 1 / sqrt x)"
      by (rule landau_o.big.compose[OF _ sqrt_at_top]) simp_all
    finally have "(λx. x * ?f x)  O(λx. x * (1 / sqrt x))" 
      by (intro landau_o.big.mult landau_o.big_refl)
    also have "(λx::real. x * (1 / sqrt x))  Θ(λx. sqrt x)"
      by (intro bigthetaI_cong eventually_mono[OF eventually_gt_at_top[of "0::real"]]) 
         (auto elim!: eventually_mono simp: field_simps)
    finally show "(λx. x * ?f x)  O(sqrt)" .
  qed fact+
  finally show ?thesis .
qed

lemma sum_upto_sqrt_bound:
  assumes x: "x  (0 :: real)"
  shows   "norm ((sum_upto (λ_. 1) (sqrt x))2 - x)  2 * norm (sqrt x)"
proof -
  from x have "0  2 * sqrt x * (1 - frac (sqrt x)) + frac (sqrt x) ^ 2"
    by (intro add_nonneg_nonneg mult_nonneg_nonneg) (simp_all add: frac_le_1)
  also from x have " = (sqrt x - frac (sqrt x)) ^ 2 - x + 2 * sqrt x"
    by (simp add: algebra_simps power2_eq_square)
  also have "sqrt x - frac (sqrt x) = of_int sqrt x" by (simp add: frac_def)
  finally have "(of_int sqrt x) ^ 2 - x  -2 * sqrt x" by (simp add: algebra_simps)
  moreover from x have "of_int (sqrt x) ^ 2  sqrt x ^ 2"
    by (intro power_mono) simp_all
  with x have "of_int (sqrt x) ^ 2 - x  0" by simp
  ultimately have "sum_upto (λ_. 1) (sqrt x) ^ 2 - x  {-2 * sqrt x..0}"
    using x by (simp add: sum_upto_altdef)
  with x show ?thesis by simp
qed

lemma summatory_divisor_count_asymptotics:
  "(λx. sum_upto (λn. real (divisor_count n)) x -
          (x * ln x + (2 * euler_mascheroni - 1) * x))  O(sqrt)"
proof -
  let ?f = "λx. x * ln x / 2 + euler_mascheroni * x"
  have "(λx. sum_upto (dirichlet_prod (λ_. 1 :: real) (λ_. 1)) x - (?f x + ?f x - x))  O(sqrt)"
    (is "?g  _")
  proof (rule hyperbola_method_bigo)
    have "eventually (λx::real. norm (sum_upto (λ_. 1) (sqrt x) ^ 2 - x)  
             2 * norm (sqrt x)) at_top"
      using eventually_ge_at_top[of "0::real"] by eventually_elim (rule sum_upto_sqrt_bound)
    thus "(λx::real. sum_upto (λ_. 1) (sqrt x) * sum_upto (λ_. 1) (sqrt x) - x)  O(sqrt)"
      by (intro bigoI[of _ 2]) (simp_all add: power2_eq_square)
  next
    show "(λx. sum_upto (λn. 1 * sum_upto (λ_. 1) (x / real n)) (sqrt x) -
                 (x * ln x / 2 + euler_mascheroni * x))  O(sqrt)"
      using divisor_count_asymptotics_aux by simp
  next
    show "(λx. sum_upto (λn. sum_upto (λ_. 1) (x / real n) * 1) (sqrt x) -
                 (x * ln x / 2 + euler_mascheroni * x))  O(sqrt)"
      using divisor_count_asymptotics_aux by simp
  qed
  also have "divisor_count n = dirichlet_prod (λ_. 1) (λ_. 1) n" for n
    using fds_divisor_count
    by (cases "n = 0") (simp_all add: fds_eq_iff power2_eq_square fds_nth_mult)
  hence "?g = (λx. sum_upto (λn. real (divisor_count n)) x - 
               (x * ln x + (2 * euler_mascheroni - 1) * x))"
    by (intro ext) (simp_all add: algebra_simps dirichlet_prod_def)
  finally show ?thesis .
qed

theorem summatory_divisor_count_asymptotics':
  "(λx. sum_upto (λn. real (divisor_count n)) x) =o 
     (λx. x * ln x + (2 * euler_mascheroni - 1) * x) +o O(λx. sqrt x)"
  using summatory_divisor_count_asymptotics
  by (subst set_minus_plus [symmetric]) (simp_all add: fun_diff_def)

theorem summatory_divisor_count_asymptotics'':
  "sum_upto (λn. real (divisor_count n)) ∼[at_top] (λx. x * ln x)"
proof -
  have "(λx. sum_upto (λn. real (divisor_count n)) x - 
           (x * ln x + (2 * euler_mascheroni - 1) * x))  O(sqrt)"
    by (rule summatory_divisor_count_asymptotics)
  also have "sqrt  Θ(λx. x powr (1/2))"
    by (intro bigthetaI_cong eventually_mono [OF eventually_ge_at_top[of "0::real"]]) 
       (auto elim!: eventually_mono simp: powr_half_sqrt)
  also have "(λx::real. x powr (1/2))  o(λx. x * ln x + (2 * euler_mascheroni - 1) * x)" by simp
  finally have "sum_upto (λn. real (divisor_count n)) ∼[at_top]
                  (λx. x * ln x + (2 * euler_mascheroni - 1) * x)"
    by (simp add: asymp_equiv_altdef)
  also have " ∼[at_top] (λx. x * ln x)" by (subst asymp_equiv_add_right) simp_all
  finally show ?thesis .
qed

lemma summatory_divisor_eq:
  "sum_upto (λn. real (divisor_count n)) (real m) = card {(n,d). n  {0<..m}  d dvd n}"
proof -
  have "sum_upto (λn. real (divisor_count n)) m = card (SIGMA n:{0<..m}. {d. d dvd n})"
    unfolding sum_upto_altdef divisor_count_def by (subst card_SigmaI) simp_all
  also have "(SIGMA n:{0<..m}. {d. d dvd n}) = {(n,d). n  {0<..m}  d dvd n}" by auto
  finally show ?thesis .
qed

context
  fixes M :: "nat  real"
  defines "M  λm. card {(n,d). n  {0<..m}  d dvd n} / card {0<..m}"
begin

lemma mean_divisor_count_asymptotics:
  "(λm. M m - (ln m + 2 * euler_mascheroni - 1))  O(λm. 1 / sqrt m)"
proof -
  have "(λm. M m - (ln m + 2 * euler_mascheroni - 1)) 
           Θ(λm. (sum_upto (λn. real (divisor_count n)) (real m) -
                 (m * ln m + (2 * euler_mascheroni - 1) * m)) / m)" (is "_  Θ(?f)")
    unfolding M_def 
    by (intro bigthetaI_cong eventually_mono [OF eventually_gt_at_top[of "0::nat"]]) 
       (auto simp: summatory_divisor_eq field_simps)
  also have "?f  O(λm. sqrt m / m)"
    by (intro landau_o.big.compose[OF _ filterlim_real_sequentially] landau_o.big.divide_right
          summatory_divisor_count_asymptotics eventually_at_top_not_equal)
  also have "(λm::nat. sqrt m / m)  Θ(λm. 1 / sqrt m)"
    by (intro bigthetaI_cong eventually_mono [OF eventually_gt_at_top[of "0::nat"]]) 
       (auto simp: field_simps)
  finally show ?thesis .
qed

theorem mean_divisor_count_asymptotics':
  "M =o (λx. ln x + 2 * euler_mascheroni - 1) +o O(λx. 1 / sqrt x)"
  using mean_divisor_count_asymptotics
  by (subst set_minus_plus [symmetric]) (simp_all add: fun_diff_def)

theorem mean_divisor_count_asymptotics'': "M ∼[at_top] ln"
proof -
  have "(λx. M x - (ln x + 2 * euler_mascheroni - 1))  O(λx. 1 / sqrt x)"
    by (rule mean_divisor_count_asymptotics)
  also have "(λx. 1 / sqrt (real x))  Θ(λx. x powr (-1/2))"
    using eventually_gt_at_top[of "0::nat"] 
    by (intro bigthetaI_cong) 
       (auto elim!: eventually_mono simp: powr_half_sqrt field_simps powr_minus)
  also have "(λx::nat. x powr (-1/2))  o(λx. ln x + 2 * euler_mascheroni - 1)"
    by (intro smallo_real_nat_transfer) simp_all
  finally have "M ∼[at_top] (λx. ln x + 2 * euler_mascheroni - 1)"
    by (simp add: asymp_equiv_altdef)
  also have " = (λx::nat. ln x + (2 * euler_mascheroni - 1))" by (simp add: algebra_simps)
  also have " ∼[at_top] (λx::nat. ln x)" by (subst asymp_equiv_add_right) auto
  finally show ?thesis .
qed

end


subsection ‹The asymptotic ditribution of coprime pairs›

context
  fixes A :: "nat  (nat × nat) set"
  defines "A  (λN. {(m,n)  {1..N} × {1..N}. coprime m n})"  
begin

lemma coprime_pairs_asymptotics:
  "(λN. real (card (A N)) - 6 / pi2 * (real N)2)  O(λN. real N * ln (real N))"
proof -
  define C :: "nat  (nat × nat) set"
    where "C = (λN. (m{1..N}. (λn. (m,n)) ` totatives m))"
  define D :: "nat  (nat × nat) set"
    where "D = (λN. (n{1..N}. (λm. (m,n)) ` totatives n))"
  have fin: "finite (C N)" "finite (D N)" for N unfolding C_def D_def
    by (intro finite_UN_I finite_imageI; simp)+    
  
  have *: "card (A N) = 2 * (m{0<..N}. totient m) - 1" if N: "N > 0" for N
  proof -
    have "A N = C N  D N"
      by (auto simp add: A_def C_def D_def totatives_def image_iff ac_simps)
    also have "card  = card (C N) + card (D N) - card (C N  D N)"
      using card_Un_Int[OF fin[of N]] by arith
    also have "C N  D N = {(1, 1)}" using N by (auto simp: image_iff totatives_def C_def D_def)
    also have "D N = (λ(x,y). (y,x)) ` C N" by (simp add: image_UN image_image C_def D_def)
    also have "card  = card (C N)" by (rule card_image) (simp add: inj_on_def C_def)
    also have "card (C N) = (m{1..N}. card ((λn. (m,n)) ` totatives m))"
      unfolding C_def by (intro card_UN_disjoint) auto
    also have " = (m{1..N}. totient m)" unfolding totient_def
      by (subst card_image) (auto simp: inj_on_def)
    also have " = (m{0<..N}. totient m)" by (intro sum.cong refl) auto
    finally show "card (A N) = 2 *  - 1" by simp
  qed
  have **: "(m{0<..N}. totient m)  1" if "N  1" for N
  proof -
    have "1  N" by fact
    also have "N = (m{0<..N}. 1)" by simp
    also have "(m{0<..N}. 1)  (m{0<..N}. totient m)" 
      by (intro sum_mono) (simp_all add: Suc_le_eq)
    finally show ?thesis .
  qed
  
  have "(λN. real (card (A N)) - 6 / pi2 * (real N)2) 
          Θ(λN. 2 * (sum_upto (λm. real (totient m)) (real N) - (3 / pi2 * (real N)2)) - 1)"
    (is "_  Θ(?f)") using * **
    by (intro bigthetaI_cong eventually_mono [OF eventually_gt_at_top[of "0::nat"]]) 
       (auto simp: of_nat_diff sum_upto_altdef)
  also have "?f  O(λN. real N * ln (real N))"
  proof (rule landau_o.big.compose[OF _ filterlim_real_sequentially], rule sum_in_bigo)
    show " (λx. 2 * (sum_upto (λm. real (totient m)) x - 3 / pi2 * x2))  O(λx. x * ln x)"
      by (subst landau_o.big.cmult_in_iff, simp, rule summatory_totient_asymptotics)
  qed simp_all
  finally show ?thesis .
qed

theorem coprime_pairs_asymptotics':
  "(λN. real (card (A N))) =o (λN. 6 / pi2 * (real N)2) +o O(λN. real N * ln (real N))"
  using coprime_pairs_asymptotics
  by (subst set_minus_plus [symmetric]) (simp_all add: fun_diff_def)

theorem coprime_pairs_asymptotics'':
  "(λN. real (card (A N))) ∼[at_top] (λN. 6 / pi2 * (real N)2)"
proof -
  have "(λN. real (card (A N)) - 6 / pi2 * (real N) ^ 2)  O(λN. real N * ln (real N))"
    by (rule coprime_pairs_asymptotics)
  also have "(λN. real N * ln (real N))  o(λN. 6 / pi ^ 2 * real N ^ 2)"
    by (rule landau_o.small.compose[OF _ filterlim_real_sequentially]) simp
  finally show ?thesis by (simp add: asymp_equiv_altdef)
qed

theorem coprime_probability_tendsto:
  "(λN. card (A N) / card ({1..N} × {1..N}))  6 / pi2"
proof -
  have "(λN. 6 / pi ^ 2) ∼[at_top] (λN. 6 / pi ^ 2 * real N ^ 2 / real N ^ 2)" 
    using eventually_gt_at_top[of "0::nat"]
    by (intro asymp_equiv_refl_ev) (auto elim!: eventually_mono)
  also have " ∼[at_top] (λN. real (card (A N)) / real N ^ 2)"
    by (intro asymp_equiv_intros asymp_equiv_symI[OF coprime_pairs_asymptotics''])
  also have " ∼[at_top] (λN. real (card (A N)) / real (card ({1..N} × {1..N})))"
    by (simp add: power2_eq_square)
  finally have " ∼[at_top] (λ_. 6 / pi ^ 2)" by (simp add: asymp_equiv_sym)
  thus ?thesis by (rule asymp_equivD_const)
qed

end


subsection ‹The asymptotics of the number of Farey fractions›

definition farey_fractions :: "nat  rat set" where
  "farey_fractions N = {q :: rat  {0<..1}. snd (quotient_of q)  int N}  "

lemma Fract_eq_coprime: 
  assumes "Rat.Fract a b = Rat.Fract c d" "b > 0" "d > 0" "coprime a b" "coprime c d"
  shows   "a = c" "b = d"
proof -
  from assms have "a * d = c * b" by (auto simp: eq_rat)
  hence "abs (a * d) = abs (c * b)" by (simp only:)
  hence "abs a * abs d = abs c * abs b" by (simp only: abs_mult)
  also have "?this  abs a = abs c  d = b" 
    using assms by (subst coprime_crossproduct_int) simp_all
  finally show "b = d" by simp
  with a * d = c * b and b > 0 show "a = c" by simp
qed
  
lemma quotient_of_split: 
  "P (quotient_of q) = (a b. b > 0  coprime a b  q = Rat.Fract a b  P (a, b))"
  by (cases q) (auto simp: quotient_of_Fract dest: Fract_eq_coprime)

lemma quotient_of_split_asm: 
  "P (Rat.quotient_of q) = (¬(a b. b > 0  coprime a b  q = Rat.Fract a b  ¬P (a, b)))"
  using quotient_of_split[of P q] by blast

lemma farey_fractions_bij:
  "bij_betw (λ(a,b). Rat.Fract (int a) (int b)) 
     {(a,b)|a b. 0 < a  a  b  b  N  coprime a b} (farey_fractions N)"
proof (rule bij_betwI[of _ _ _ "λq. case quotient_of q of (a, b)  (nat a, nat b)"], goal_cases)
  case 1
  show ?case
    by (auto simp: farey_fractions_def Rat.zero_less_Fract_iff Rat.Fract_le_one_iff 
                   Rat.quotient_of_Fract Rat.normalize_def gcd_int_def Let_def)
next
  case 2
  show ?case
    by (auto simp add: farey_fractions_def Rat.Fract_le_one_iff Rat.zero_less_Fract_iff split: prod.splits quotient_of_split_asm)
      (simp add: coprime_int_iff [symmetric])
next
  case (3 x)
  thus ?case by (auto simp: Rat.quotient_of_Fract Rat.normalize_def Let_def gcd_int_def)
next
  case (4 x)
  thus ?case unfolding farey_fractions_def
    by (split quotient_of_split) (auto simp: Rat.zero_less_Fract_iff)
qed

lemma card_farey_fractions: "card (farey_fractions N) = sum totient {0<..N}"
proof -
  have "card (farey_fractions N) = card {(a,b)|a b. 0 < a  a  b  b  N  coprime a b}"
    using farey_fractions_bij by (rule bij_betw_same_card [symmetric])
  also have "{(a,b)|a b. 0 < a  a  b  b  N  coprime a b} =
               (b{0<..N}. (λa. (a, b)) ` totatives b)"
    by (auto simp: totatives_def image_iff)
  also have "card  = (b{0<..N}. card ((λa. (a, b)) ` totatives b))"
    by (intro card_UN_disjoint) auto
  also have " = (b{0<..N}. totient b)"
    unfolding totient_def by (intro sum.cong refl card_image) (auto simp: inj_on_def)
  finally show ?thesis .
qed

lemma card_farey_fractions_asymptotics:
  "(λN. real (card (farey_fractions N)) - 3 / pi2 * (real N)2)  O(λN. real N * ln (real N))"
proof -
  have "(λN. sum_upto (λn. real (totient n)) (real N) - 3 / pi2 * (real N)2) 
             O(λN. real N * ln (real N))" (is "?f  _")
    using summatory_totient_asymptotics filterlim_real_sequentially 
    by (rule landau_o.big.compose)
  also have "?f = (λN. real (card (farey_fractions N)) - 3 / pi2 * (real N)2)"
    by (intro ext) (simp add: sum_upto_altdef card_farey_fractions)
  finally show ?thesis .
qed

theorem card_farey_fractions_asymptotics':
  "(λN. card (farey_fractions N)) =o (λN. 3 / pi2 * N^2) +o O(λN. N * ln N)"
  using card_farey_fractions_asymptotics
  by (subst set_minus_plus [symmetric]) (simp_all add: fun_diff_def)

theorem card_farey_fractions_asymptotics'':
  "(λN. real (card (farey_fractions N))) ∼[at_top] (λN. 3 / pi2 * (real N)2)"
proof -
  have "(λN. real (card (farey_fractions N)) - 3 / pi2 * (real N) ^ 2)  O(λN. real N * ln (real N))"
    by (rule card_farey_fractions_asymptotics)
  also have "(λN. real N * ln (real N))  o(λN. 3 / pi ^ 2 * real N ^ 2)"
    by (rule landau_o.small.compose[OF _ filterlim_real_sequentially]) simp
  finally show ?thesis by (simp add: asymp_equiv_altdef)
qed

end