Session Prime_Harmonic_Series

Theory Prime_Harmonic_Misc

  File:   Prime_Harmonic_Misc.thy
  Author: Manuel Eberl <>


section ‹Auxiliary lemmas›
theory Prime_Harmonic_Misc

lemma sum_list_nonneg: "xset xs. x  0  sum_list xs  (0 :: 'a :: ordered_ab_group_add)"
  by (induction xs) auto

lemma sum_telescope':
  assumes "m  n"
  shows   "(k = Suc m..n. f k - f (Suc k)) = f (Suc m) - (f (Suc n) :: 'a :: ab_group_add)"
  by (rule dec_induct[OF assms]) (simp_all add: algebra_simps)

lemma dvd_prodI:
  assumes "finite A" "x  A"
  shows   "f x dvd prod f A"
proof -
  from assms have "prod f A = f x * prod f (A - {x})" 
    by (intro prod.remove) simp_all
  thus ?thesis by simp

lemma dvd_prodD: "finite A  prod f A dvd x  a  A  f a dvd x"
  by (erule dvd_trans[OF dvd_prodI])

lemma multiplicity_power_nat: 
  "prime p  n > 0  multiplicity p (n ^ k :: nat) = k * multiplicity p n"
  by (induction k) (simp_all add: prime_elem_multiplicity_mult_distrib)

lemma multiplicity_prod_prime_powers_nat':
  "finite S  pS. prime p  prime p  
     multiplicity p (S :: nat) = (if p  S then 1 else 0)"
  using multiplicity_prod_prime_powers[of S p "λ_. 1"] by simp

lemma prod_prime_subset:
  assumes "finite A" "finite B"
  assumes "x. x  A  prime (x::nat)"
  assumes "x. x  B  prime x"
  assumes "A dvd B"
  shows   "A  B"
  fix x assume x: "x  A"
  from assms(4)[of 0] have "0  B" by auto
  with assms have nonzero: "zB. z  0" by (intro ballI notI) auto

  from x assms have "1 = multiplicity x (A)"
    by (subst multiplicity_prod_prime_powers_nat') simp_all
  also from assms nonzero have "  multiplicity x (B)" by (intro dvd_imp_multiplicity_le) auto
  finally have "multiplicity x (B) > 0" by simp
  moreover from assms x have "prime x" by simp
  ultimately show "x  B" using assms(2,4)
    by (subst (asm) multiplicity_prod_prime_powers_nat') (simp_all split: if_split_asm)

lemma prod_prime_eq:
  assumes "finite A" "finite B" "x. x  A  prime (x::nat)" "x. x  B  prime x" "A = B"
  shows   "A = B"
  using assms by (intro equalityI prod_prime_subset) simp_all

lemma ln_ln_nonneg:
  assumes x: "x  (3 :: real)"
  shows   "ln (ln x)  0"
proof -
  have "exp 1  (3::real)" by (rule  exp_le)
  hence "ln (exp 1)  ln (3 :: real)" by (subst ln_le_cancel_iff) simp_all
  also from x have "  ln x" by (subst ln_le_cancel_iff) simp_all
  finally have "ln 1  ln (ln x)" using x by (subst ln_le_cancel_iff) simp_all
  thus ?thesis by simp


Theory Squarefree_Nat

  File:   Squarefree_Nat.thy
  Author: Manuel Eberl <>

  The unique decomposition of a natural number into a squarefree part and a square.

section ‹Squarefree decomposition of natural numbers›
theory Squarefree_Nat

text ‹
  The squarefree part of a natural number is the set of all prime factors that appear 
  with odd multiplicity. The square part, correspondingly, is what remains after dividing
  by the squarefree part.
definition squarefree_part :: "nat  nat set" where
  "squarefree_part n = {pprime_factors n. odd (multiplicity p n)}"

definition square_part :: "nat  nat" where
  "square_part n = (if n = 0 then 0 else (pprime_factors n. p ^ (multiplicity p n div 2)))"

lemma squarefree_part_0 [simp]: "squarefree_part 0 = {}"
  by (simp add: squarefree_part_def)

lemma square_part_0 [simp]: "square_part 0 = 0"
  by (simp add: square_part_def)

lemma squarefree_decompose: "(squarefree_part n) * square_part n ^ 2 = n"
proof (cases "n = 0")
  case False
  define A s where "A = squarefree_part n" and "s = square_part n"
  have "(A) = (pA. p ^ (multiplicity p n mod 2))"
    by (intro prod.cong) (auto simp: A_def squarefree_part_def elim!: oddE)
  also have " = (pprime_factors n. p ^ (multiplicity p n mod 2))"
    by (intro prod.mono_neutral_left) (auto simp: A_def  squarefree_part_def)
  also from False have " * s^2 = n"
    by (simp add: s_def square_part_def prod.distrib [symmetric] power_add [symmetric] 
                  power_mult [symmetric] prime_factorization_nat [symmetric] algebra_simps
  finally show "A * s^2 = n" .
qed simp

lemma squarefree_part_pos [simp]: "(squarefree_part n) > 0"
  using prime_gt_0_nat unfolding squarefree_part_def by auto
lemma squarefree_part_ge_Suc_0 [simp]: "(squarefree_part n)  Suc 0"
  using squarefree_part_pos[of n] by presburger

lemma squarefree_part_subset [intro]: "squarefree_part n  prime_factors n"
  unfolding squarefree_part_def by auto

lemma squarefree_part_finite [simp]: "finite (squarefree_part n)"
  by (rule finite_subset[OF squarefree_part_subset]) simp

lemma squarefree_part_dvd [simp]: "(squarefree_part n) dvd n"
  by (subst (2) squarefree_decompose [of n, symmetric]) simp

lemma squarefree_part_dvd' [simp]: "p  squarefree_part n  p dvd n"
  by (rule dvd_prodD[OF _ squarefree_part_dvd]) simp_all

lemma square_part_dvd [simp]: "square_part n ^ 2 dvd n"
  by (subst (2) squarefree_decompose [of n, symmetric]) simp

lemma square_part_dvd' [simp]: "square_part n dvd n"
  by (subst (2) squarefree_decompose [of n, symmetric]) simp

lemma squarefree_part_le: "p  squarefree_part n  p  n"
  by (cases "n = 0") (simp_all add: dvd_imp_le)

lemma square_part_le: "square_part n  n"
  by (cases "n = 0") (simp_all add: dvd_imp_le)

lemma square_part_le_sqrt: "square_part n  nat sqrt (real n)"
proof -
  have "1 * square_part n ^ 2  (squarefree_part n) * square_part n ^ 2"
    by (intro mult_right_mono) simp_all
  also have " = n" by (rule squarefree_decompose)
  finally have "real (square_part n ^ 2)  real n" by (subst of_nat_le_iff) simp
  hence "sqrt (real (square_part n ^ 2))  sqrt (real n)" 
    by (subst real_sqrt_le_iff) simp_all
  also have "sqrt (real (square_part n ^ 2)) = real (square_part n)" by simp
  finally show ?thesis by linarith

lemma square_part_pos [simp]: "n > 0  square_part n > 0"
  unfolding square_part_def using prime_gt_0_nat by auto
lemma square_part_ge_Suc_0 [simp]: "n > 0  square_part n  Suc 0"
  using square_part_pos[of n] by presburger

lemma zero_not_in_squarefree_part [simp]: "0  squarefree_part n"
  unfolding squarefree_part_def by auto

lemma multiplicity_squarefree_part:
  "prime p  multiplicity p ((squarefree_part n)) = (if p  squarefree_part n then 1 else 0)"
  using squarefree_part_subset[of n]
  by (intro multiplicity_prod_prime_powers_nat') auto

text ‹
  The squarefree part really is square, its only square divisor is 1.
lemma square_dvd_squarefree_part_iff:
  "x^2 dvd (squarefree_part n)  x = 1"
proof (rule iffI, rule multiplicity_eq_nat)
  assume dvd: "x^2 dvd (squarefree_part n)"
  hence "x  0" using squarefree_part_subset[of n] prime_gt_0_nat by (intro notI) auto
  thus x: "x > 0" by simp
  fix p :: nat assume p: "prime p"
  from p x have "2 * multiplicity p x = multiplicity p (x^2)" 
    by (simp add: multiplicity_power_nat)
  also from dvd have "  multiplicity p ((squarefree_part n))"
    by (intro dvd_imp_multiplicity_le) simp_all
  also have "  1" using multiplicity_squarefree_part[of p n] p by simp
  finally show "multiplicity p x = multiplicity p 1" by simp
qed simp_all

lemma squarefree_decomposition_unique1:
  assumes "squarefree_part m = squarefree_part n"
  assumes "square_part m = square_part n"
  shows   "m = n"
  by (subst (1 2) squarefree_decompose [symmetric]) (simp add: assms)

lemma squarefree_decomposition_unique2:
  assumes n: "n > 0"
  assumes decomp: "n = (A2 * s2^2)"
  assumes prime: "x. x  A2  prime x"
  assumes fin: "finite A2"
  assumes s2_nonneg: "s2  0"
  shows "A2 = squarefree_part n" and "s2 = square_part n"
proof -
  define A1 s1 where "A1 = squarefree_part n" and "s1 = square_part n"
  have "finite A1" unfolding A1_def by simp
  note fin = ‹finite A1 ‹finite A2

  have "A1  prime_factors n" unfolding A1_def using squarefree_part_subset .
  note subset = this prime

  have "A1 > 0" "A2 > 0" using subset(1)  prime_gt_0_nat 
    by (auto intro!: prod_pos dest: prime)
  from n have "s1 > 0" unfolding s1_def by simp
  from assms have "s2  0" by (intro notI) simp
  hence "s2 > 0" by simp
  note pos = A1 > 0 A2 > 0 s1 > 0 s2 > 0

  have eq': "multiplicity p s1 = multiplicity p s2" 
            "multiplicity p (A1) = multiplicity p (A2)" 
    if   p: "prime p" for p
  proof -
    define m where "m = multiplicity p"
    from decomp have "m (A1 * s1^2) = m (A2 * s2^2)" unfolding A1_def s1_def
      by (simp add: A1_def s1_def squarefree_decompose)
    with p pos have eq: "m (A1) + 2 * m s1 = m (A2) + 2 * m s2"
      by (simp add: m_def prime_elem_multiplicity_mult_distrib multiplicity_power_nat)
    moreover from fin subset p have "m (A1)  1" "m (A2)  1" unfolding m_def
      by ((subst multiplicity_prod_prime_powers_nat', auto)[])+
    ultimately show "m s1 = m s2" by linarith
    with eq show "m (A1) = m (A2)" by simp
  show "s2 = square_part n"
    by (rule multiplicity_eq_nat) (insert pos eq'(1), auto simp: s1_def)
  have "A2 = (squarefree_part n)"
    by (rule multiplicity_eq_nat) (insert pos eq'(2), auto simp: A1_def)
  with fin subset show "A2 = squarefree_part n" unfolding A1_def
    by (intro prod_prime_eq) auto

lemma squarefree_decomposition_unique2':
  assumes decomp: "(A1 * s1^2) = (A2 * s2^2 :: nat)"
  assumes fin: "finite A1" "finite A2"
  assumes subset: "x. x  A1  prime x" "x. x  A2  prime x"
  assumes pos: "s1 > 0" "s2 > 0"
  defines "n  A1 * s1^2"
  shows "A1 = A2" "s1 = s2"
proof -
  from pos have n: "n > 0" using prime_gt_0_nat 
    by (auto simp: n_def intro!: prod_pos dest: subset)
  have "A1 = squarefree_part n" "s1 = square_part n" 
    by ((rule squarefree_decomposition_unique2[of n A1 s1], insert assms n, simp_all)[])+
  moreover have "A2 = squarefree_part n" "s2 = square_part n" 
    by ((rule squarefree_decomposition_unique2[of n A2 s2], insert assms n, simp_all)[])+
  ultimately show "A1 = A2" "s1 = s2" by simp_all

text ‹
  The following is a nice and simple lower bound on the number of prime numbers less than 
  a given number due to Erd\H{o}s. In particular, it implies that there are infinitely many primes.
lemma primes_lower_bound:
  fixes n :: nat
  assumes "n > 0"
  defines "π  λn. card {p. prime p  p  n}"
  shows   "real (π n)  ln (real n) / ln 4"
proof -
  have "real n = real (card {1..n})" by simp
  also have "{1..n} = (λ(A, b). A * b^2) ` (λn. (squarefree_part n, square_part n)) ` {1..n}"
    unfolding image_comp o_def squarefree_decompose case_prod_unfold fst_conv snd_conv by simp
  also have "card   card ((λn. (squarefree_part n, square_part n)) ` {1..n})"
    by (rule card_image_le) simp_all
  also have "  card (squarefree_part ` {1..n} × square_part ` {1..n})"
    by (rule card_mono) auto
  also have "real  = real (card (squarefree_part ` {1..n})) * real (card (square_part ` {1..n}))"
    by simp
  also have "  2 ^ π n * sqrt (real n)"
  proof (rule mult_mono)
    have "card (squarefree_part ` {1..n})  card (Pow {p. prime p  p  n})"
      using squarefree_part_subset squarefree_part_le by (intro card_mono) force+
    also have "real  = 2 ^ π n" by (simp add: π_def card_Pow)
    finally show "real (card (squarefree_part ` {1..n}))  2 ^ π n" by - simp_all
    have "square_part k  nat sqrt n" if "k  n" for k
      by (rule order.trans[OF square_part_le_sqrt])
         (insert that, auto intro!: nat_mono floor_mono)
    hence "card (square_part ` {1..n})  card {1..nat sqrt n}"
      by (intro card_mono) (auto intro: order.trans[OF square_part_le_sqrt])
    also have " = nat sqrt n" by simp
    also have "real   sqrt n" by simp
    finally show "real (card (square_part ` {1..n}))  sqrt (real n)" by - simp_all
  qed simp_all
  finally have "real n  2 ^ π n * sqrt (real n)" by - simp_all
  with n > 0 have "ln (real n)  ln (2 ^ π n * sqrt (real n))"
    by (subst ln_le_cancel_iff) simp_all
  moreover have "ln (4 :: real) = real 2 * ln 2" by (subst ln_realpow [symmetric]) simp_all
  ultimately show ?thesis using n > 0
    by (simp add: ln_mult ln_realpow[of _ "π n"] ln_sqrt field_simps)


Theory Prime_Harmonic

  File:   Prime_Harmonic.thy
  Author: Manuel Eberl <>

  A lower bound for the partial sums of the prime harmonic series, and a proof of its divergence.
  (#81 on the list of 100 mathematical theorems)

section ‹The Prime Harmonic Series›
theory Prime_Harmonic

subsection ‹Auxiliary equalities and inequalities›

text ‹
  First of all, we prove the following result about rearranging a product over a set into a sum
  over all subsets of that set.
lemma prime_harmonic_aux1:
  fixes A :: "'a :: field set"
  shows "finite A  (xA. 1 + 1 / x) = (xPow A. 1 / x)"
proof (induction rule: finite_induct)
  fix a :: 'a and A :: "'a set"
  assume a: "a  A" and fin: "finite A"
  assume IH: "(xA. 1 + 1 / x) = (xPow A. 1 / x)"
  from a and fin have "(xinsert a A. 1 + 1 / x) = (1 + 1 / a) * (xA. 1 + 1 / x)" by simp
  also from fin have " = (xPow A. 1 / x) + (xPow A. 1 / (a * x))"
    by (subst IH) (auto simp add: algebra_simps sum_divide_distrib)
  also from fin a have "(xPow A. 1 / (a * x)) = (xPow A. 1 / (insert a x))"
    by (intro sum.cong refl, subst prod.insert) (auto dest: finite_subset)
  also from a have " = (xinsert a ` Pow A. 1 / x)"
    by (subst sum.reindex) (auto simp: inj_on_def)
  also from fin a have "(xPow A. 1 / x) +  = (xPow A  insert a ` Pow A. 1 / x)"
    by (intro sum.union_disjoint [symmetric]) (simp, simp, blast)
  also have "Pow A  insert a ` Pow A = Pow (insert a A)" by (simp only: Pow_insert)
  finally show " (xinsert a A. 1 + 1 / x) = (xPow (insert a A). 1 / x)" .
qed simp

text ‹
  Next, we prove a simple and reasonably accurate upper bound for the sum of the squares of any
  subset of the natural numbers, derived by simple telescoping. Our upper bound is approximately
  1.67; the exact value is $\frac{\pi^2}{6} \approx 1.64$. (cf. Basel problem)
lemma prime_harmonic_aux2:
  assumes "finite (A :: nat set)"
  shows   "(kA. 1 / (real k ^ 2))  5/3"
proof -
  define n where "n = max 2 (Max A)"
  have n: "n  Max A" "n  2" by (auto simp: n_def)
  with assms have "A  {0..n}" by (auto intro: order.trans[OF Max_ge])
  hence "(kA. 1 / (real k ^ 2))  (k=0..n. 1 / (real k ^ 2))" by (intro sum_mono2) auto
  also from n have " = 1 + (k=Suc 1..n. 1 / (real k ^ 2))" by (simp add: sum.atLeast_Suc_atMost)
  also have "(k=Suc 1..n. 1 / (real k ^ 2)) 
          (k=Suc 1..n. 1 / (real k ^ 2 - 1/4))" unfolding power2_eq_square
    by (intro sum_mono divide_left_mono mult_pos_pos)
       (linarith, simp_all add: field_simps less_1_mult)
  also have " = (k=Suc 1..n. 1 / (real k - 1/2) - 1 / (real (Suc k) - 1/2))"
    by (intro sum.cong refl) (simp_all add: field_simps power2_eq_square)
  also from n have " = 2 / 3 - 1 / (1 / 2 + real n)"
    by (subst sum_telescope') simp_all
  also have "1 +   5/3" by simp
  finally show ?thesis by - simp

subsection ‹Estimating the partial sums of the Prime Harmonic Series›

text ‹
  We are now ready to show our main result: the value of the partial prime harmonic sum over
  all primes no greater than $n$ is bounded from below by the $n$-th harmonic number
  $H_n$ minus some constant.

  In our case, this constant will be $\frac{5}{3}$. As mentioned before, using a
  proof of the Basel problem can improve this to $\frac{\pi^2}{6}$, but the improvement is very
  small and the proof of the Basel problem is a very complex one.

  The exact asymptotic behaviour of the partial sums is actually $\ln (\ln n) + M$, where $M$
  is the Meissel--Mertens constant (approximately 0.261).
theorem prime_harmonic_lower:
  assumes n: "n  2"
  shows "(pprimes_upto n. 1 / real p)  ln (harm n) - ln (5/3)"
proof -
  ― ‹the set of primes that we will allow in the squarefree part›
  define P where "P n = set (primes_upto n)" for n
    fix n :: nat
    have "finite (P n)" by (simp add: P_def)
  } note [simp] = this

  ― ‹The function that combines the squarefree part and the square part›
  define f where "f = (λ(R, s :: nat). R * s^2)"

  ― ‹@{term f} is injective if the squarefree part contains only primes
      and the square part is positive.›
  have inj: "inj_on f (Pow (P n)×{1..n})"
  proof (rule inj_onI, clarify, rule conjI)
    fix A1 A2 :: "nat set" and s1 s2 :: nat
    assume A: "A1  P n" "A2  P n" "s1  {1..n}" "s2  {1..n}" "f (A1, s1) = f (A2, s2)"
    have fin: "finite A1" "finite A2" by (rule A(1,2)[THEN finite_subset], simp)+
    show "A1 = A2" "s1 = s2"
      by ((rule squarefree_decomposition_unique2'[of A1 s1 A2 s2],
          insert A fin, auto simp: f_def P_def set_primes_upto)[])+

  ― ‹@{term f} hits every number between @{term "1::nat"} and @{term "n"}. It also hits a lot
      of other numbers, but we do not care about those, since we only need a lower bound.›
  have surj: "{1..n}  f ` (Pow (P n)×{1..n})"
    fix x assume x: "x  {1..n}"
    have "x = f (squarefree_part x, square_part x)" by (simp add: f_def squarefree_decompose)
    moreover have "squarefree_part x  Pow (P n)" using squarefree_part_subset[of x] x
      by (auto simp: P_def set_primes_upto intro: order.trans[OF squarefree_part_le[of _ x]])
    moreover have "square_part x  {1..n}" using x
      by (auto simp: Suc_le_eq intro: order.trans[OF square_part_le[of x]])
    ultimately show "x  f ` (Pow (P n)×{1..n})" by simp

  ― ‹We now show the main result by rearranging the sum over all primes to a product over all
      all squarefree parts times a sum over all square parts, and then applying some simple-minded
  have "harm n = (n=1..n. 1 / real n)" by (simp add: harm_def field_simps)
  also from surj have "  (nf ` (Pow (P n)×{1..n}). 1 / real n)"
    by (intro sum_mono2 finite_imageI finite_cartesian_product) simp_all
  also from inj have " = (xPow (P n)×{1..n}. 1 / real (f x))"
    by (subst sum.reindex) simp_all
  also have " = (APow (P n). 1 / real (A)) * (k=1..n. 1 / (real k)^2)" unfolding f_def
    by (subst sum_product, subst sum.cartesian_product) (simp add: case_prod_beta)
  also have "  (APow (P n). 1 / real (A)) * (5/3)"
    by (intro mult_left_mono prime_harmonic_aux2 sum_nonneg)
       (auto simp: P_def intro!: prod_nonneg)
  also have "(APow (P n). 1 / real (A)) = (A((`) real) ` Pow (P n). 1 / A)"
    by (subst sum.reindex) (auto simp: inj_on_def inj_image_eq_iff prod.reindex)
  also have "((`) real) ` Pow (P n) = Pow (real ` P n)" by (intro image_Pow_surj refl)
  also have "(APow (real ` P n). 1 / A) = (xreal ` P n. 1 + 1 / x)"
    by (intro prime_harmonic_aux1 [symmetric] finite_imageI) simp_all
  also have " = (iP n. 1 + 1 / real i)" by (subst prod.reindex) (auto simp: inj_on_def)
  also have "  (iP n. exp (1 / real i))" by (intro prod_mono) auto
  also have " = exp (iP n. 1 / real i)" by (simp add: exp_sum)
  finally have "ln (harm n)  ln ( * (5/3))" using n
    by (subst ln_le_cancel_iff) simp_all
  hence "ln (harm n) - ln (5/3)  (iP n. 1 / real i)"
    by (subst (asm) ln_mult) (simp_all add: algebra_simps)
  thus ?thesis unfolding P_def
    by (subst (asm) sum.distinct_set_conv_list) simp_all

text ‹
  We can use the inequality $\ln (n + 1) \le H_n$ to estimate the asymptotic growth of the partial
  prime harmonic series. Note that $H_n \sim \ln n + \gamma$ where $\gamma$ is the
  Euler--Mascheroni constant (approximately 0.577), so we lose some accuracy here.
corollary prime_harmonic_lower':
  assumes n: "n  2"
  shows "(pprimes_upto n. 1 / real p)  ln (ln (n + 1)) - ln (5/3)"
proof -
  from assms ln_le_harm[of n] have "ln (ln (real n + 1))  ln (harm n)" by simp
  also from assms have " - ln (5/3)  (pprimes_upto n. 1 / real p)"
    by (rule prime_harmonic_lower)
  finally show ?thesis by - simp

(* TODO: Not needed in Isabelle 2016 *)
lemma Bseq_eventually_mono:
  assumes "eventually (λn. norm (f n)  norm (g n)) sequentially" "Bseq g"
  shows   "Bseq f"
proof -
  from assms(1) obtain N where N: "n. n  N  norm (f n)  norm (g n)"
    by (auto simp: eventually_at_top_linorder)
  from assms(2) obtain K where K: "n. norm (g n)  K" by (blast elim!: BseqE)
    fix n :: nat
    have "norm (f n)  max K (Max {norm (f n) |n. n < N})"
      apply (cases "n < N")
      apply (rule max.coboundedI2, rule Max.coboundedI, auto) []
      apply (rule max.coboundedI1, force intro: order.trans[OF N K])
  thus ?thesis by (blast intro: BseqI')

lemma Bseq_add:
  assumes "Bseq (f :: nat  'a :: real_normed_vector)"
  shows   "Bseq (λx. f x + c)"
proof -
  from assms obtain K where K: "x. norm (f x)  K" unfolding Bseq_def by blast
    fix x :: nat
    have "norm (f x + c)  norm (f x) + norm c" by (rule norm_triangle_ineq)
    also have "norm (f x)  K" by (rule K)
    finally have "norm (f x + c)  K + norm c" by simp
  thus ?thesis by (rule BseqI')

lemma convergent_imp_Bseq: "convergent f  Bseq f"
  by (simp add: Cauchy_Bseq convergent_Cauchy)

(* END TODO *)

text ‹
  We now use our last estimate to show that the prime harmonic series diverges. This is obvious,
  since it is bounded from below by $\ln (\ln (n + 1))$ minus some constant, which obviously
  tends to infinite.

  Directly using the divergence of the harmonic series would also be possible and shorten this
  proof a bit..
corollary prime_harmonic_series_unbounded:
  "¬Bseq (λn. pprimes_upto n. 1 / p)" (is "¬Bseq ?f")
  assume "Bseq ?f"
  hence "Bseq (λn. ?f n + ln (5/3))" by (rule Bseq_add)
  have "Bseq (λn. ln (ln (n + 1)))"
  proof (rule Bseq_eventually_mono)
    from eventually_ge_at_top[of "2::nat"]
      show "eventually (λn. norm (ln (ln (n + 1)))  norm (?f n + ln (5/3))) sequentially"
    proof eventually_elim
      fix n :: nat assume n: "n  2"
      hence "norm (ln (ln (real n + 1))) = ln (ln (real n + 1))"
        using ln_ln_nonneg[of "real n + 1"] by simp
      also have "  ?f n + ln (5/3)" using prime_harmonic_lower'[OF n]
        by (simp add: algebra_simps)
      also have "?f n + ln (5/3)  0" by (intro add_nonneg_nonneg sum_list_nonneg) simp_all
      hence "?f n + ln (5/3) = norm (?f n + ln (5/3))" by simp
      finally show "norm (ln (ln (n + 1)))  norm (?f n + ln (5/3))"
        by (simp add: add_ac)
  qed fact
  then obtain k where k: "k > 0" "n. norm (ln (ln (real (n::nat) + 1)))  k"
    by (auto elim!: BseqE simp: add_ac)

  define N where "N = nat exp (exp k)"
  have N_pos: "N > 0" unfolding N_def by simp
  have "real N + 1 > exp (exp k)" unfolding N_def by linarith
  hence "ln (real N + 1) > ln (exp (exp k))" by (subst ln_less_cancel_iff) simp_all
  with N_pos have "ln (ln (real N + 1)) > ln (exp k)" by (subst ln_less_cancel_iff) simp_all
  hence "k < ln (ln (real N + 1))" by simp
  also have "  norm (ln (ln (real N + 1)))" by simp
  finally show False using k(2)[of N] by simp

corollary prime_harmonic_series_diverges:
  "¬convergent (λn. pprimes_upto n. 1 / p)"
  using prime_harmonic_series_unbounded convergent_imp_Bseq by blast