import Mathlib

/-!
# Divergence of the harmonic series

This file contains a self-contained formalization of the classical fact that the harmonic
series `∑ 1/n` diverges, in three equivalent forms:

* `harmonic_series_not_summable` : the family `n ↦ 1/n` is not summable over `ℕ`;
* `harmonic_partial_sums_tendsto_atTop` : the partial sums `∑_{i < N} 1/(i+1)` tend to `+∞`;
* `harmonic_partial_sums_unbounded` : the partial sums exceed every real bound.

In addition, quantifying the (very slow) rate of divergence, we prove

* `harmonic_hundred_needs_more_than_ten_pow_43_terms` : if the `n`-th harmonic number
  `H n = ∑_{i < n} 1/(i+1)` is at least `100`, then `n > 10 ^ 43`.

The last statement uses the bound `H n ≤ log n + (3/2 - log 2)` valid for `n ≥ 2`, which
follows from the monotonicity of the sequence `n ↦ H n - log n` (`H n - log n` decreases to the
Euler–Mascheroni constant `γ ≈ 0.5772`).
-/

open Filter Finset

/-- **The harmonic series diverges**: the family `n ↦ 1 / n` is not summable over `ℕ`. -/
theorem harmonic_series_not_summable : ¬ Summable (fun n : ℕ => (1 : ℝ) / n) := by
  simpa [one_div] using Real.not_summable_natCast_inv

/-- **The harmonic series diverges**: its partial sums tend to `+∞`. -/
theorem harmonic_partial_sums_tendsto_atTop :
    Filter.Tendsto (fun N : ℕ => ∑ i ∈ Finset.range N, (1 : ℝ) / (i + 1)) atTop atTop :=
  Real.tendsto_sum_range_one_div_nat_succ_atTop

/-- **The harmonic series diverges**: its partial sums exceed every real number. -/
theorem harmonic_partial_sums_unbounded (C : ℝ) :
    ∃ N : ℕ, C < ∑ i ∈ Finset.range N, (1 : ℝ) / (i + 1) := by
  obtain ⟨N, hN⟩ := (harmonic_partial_sums_tendsto_atTop.eventually_gt_atTop C).exists
  exact ⟨N, hN⟩

section Quantitative

/-- A crude but sufficient numerical lower bound for `exp 98.5`. -/
private lemma ten_pow_43_lt_two_mul_exp : (10 : ℝ) ^ 43 < 2 * Real.exp 98.5 := by
  have hpos : (0 : ℝ) < 2 * Real.exp 98.5 := by positivity
  have he : (2.7182818283 : ℝ) ^ 197 < Real.exp 197 := by
    have h1 : (2.7182818283 : ℝ) < Real.exp 1 := Real.exp_one_gt_d9
    calc (2.7182818283 : ℝ) ^ 197 < (Real.exp 1) ^ 197 := by
          exact pow_lt_pow_left₀ h1 (by norm_num) (by norm_num)
      _ = Real.exp 197 := by
          rw [← Real.exp_nat_mul]; norm_num
  have hsq : ((10 : ℝ) ^ 43) ^ 2 < (2 * Real.exp 98.5) ^ 2 := by
    have hexp : (Real.exp 98.5) ^ 2 = Real.exp 197 := by
      rw [sq, ← Real.exp_add]; norm_num
    have hnum : (10 : ℝ) ^ 86 < 4 * (2.7182818283 : ℝ) ^ 197 := by norm_num
    calc ((10 : ℝ) ^ 43) ^ 2 = (10 : ℝ) ^ 86 := by ring
      _ < 4 * (2.7182818283 : ℝ) ^ 197 := hnum
      _ < 4 * Real.exp 197 := by linarith
      _ = (2 * Real.exp 98.5) ^ 2 := by rw [mul_pow, hexp]; norm_num
  exact lt_of_pow_lt_pow_left₀ 2 hpos.le hsq

/-- For `n ≥ 2`, the `n`-th harmonic number satisfies `H n ≤ log n + (3/2 - log 2)`. -/
private lemma harmonic_le_log_add (n : ℕ) (hn : 2 ≤ n) :
    (harmonic n : ℝ) ≤ Real.log n + (3 / 2 - Real.log 2) := by
  have h := Real.strictAnti_eulerMascheroniSeq'.antitone hn
  rw [Real.eulerMascheroniSeq', Real.eulerMascheroniSeq'] at h
  have hn0 : n ≠ 0 := by omega
  simp only [hn0, if_false, show (2 : ℕ) ≠ 0 by norm_num] at h
  have h2 : (harmonic 2 : ℝ) = 3 / 2 := by norm_num [harmonic]
  rw [h2] at h
  push_cast at h ⊢
  linarith

/-- **Quantitative slowness of the divergence.** If the `n`-th harmonic number
`H n = ∑_{i < n} 1/(i+1)` reaches `100`, then `n > 10 ^ 43`. -/
theorem harmonic_hundred_needs_more_than_ten_pow_43_terms (n : ℕ) (hn : 100 ≤ harmonic n) :
    10 ^ 43 < n := by
  -- First, `n ≥ 2`.
  have hn2 : 2 ≤ n := by
    by_contra h
    interval_cases n <;> norm_num [harmonic] at hn
  have hnpos : (0 : ℝ) < n := by positivity
  have hR : (100 : ℝ) ≤ (harmonic n : ℝ) := by exact_mod_cast hn
  have hbound := harmonic_le_log_add n hn2
  -- Hence `log n ≥ 98.5 + log 2`.
  have hlog : (98.5 : ℝ) + Real.log 2 ≤ Real.log n := by linarith
  have hexp : Real.exp (98.5 + Real.log 2) ≤ (n : ℝ) := by
    calc Real.exp (98.5 + Real.log 2) ≤ Real.exp (Real.log n) := Real.exp_le_exp.2 hlog
      _ = (n : ℝ) := Real.exp_log hnpos
  have heq : Real.exp ((98.5 : ℝ) + Real.log 2) = 2 * Real.exp 98.5 := by
    rw [Real.exp_add, Real.exp_log (by norm_num : (0:ℝ) < 2), mul_comm]
  rw [heq] at hexp
  have : (10 : ℝ) ^ 43 < (n : ℝ) := lt_of_lt_of_le ten_pow_43_lt_two_mul_exp hexp
  exact_mod_cast this

end Quantitative