1. 1
2. 2
3. 3
4. 4
5. 5
6. 6
7. 7
8. 8
9. 9
10. 10
11. 11
12. 12
13. 13
14. 14
15. 15
16. 16
17. 17
18. 18
19. 19
20. 20
21. 21
22. 22
23. 23
24. 24
25. 25
26. 26
27. 27
28. 28
29. 29
30. 30
31. 31
32. 32
33. 33
34. 34
35. 35
36. 36
37. 37
38. 38
39. 39
40. 40
41. 41
42. 42
43. 43
44. 44
45. 45
46. 46
47. 47
48. 48
49. 49
50. 50
51. 51
52. 52
53. 53
54. 54
55. 55
56. 56
57. 57
58. 58
59. 59
60. 60
61. 61
62. 62
63. 63
64. 64
65. 65
66. 66
67. 67
68. 68
69. 69
70. 70
71. 71
72. 72
73. 73
74. 74
75. 75
76. 76
77. 77
78. 78
79. 79
80. 80
81. 81
82. 82
83. 83
84. 84
85. 85
86. 86
87. 87
88. 88
89. 89
90. 90
91. 91
92. 92
93. 93
94. 94
95. 95
96. 96
97. 97
98. 98
99. 99
100. 100
101. 101

import Mathlib.Algebra.Field.GeomSum
import Mathlib.Analysis.RCLike.Basic
import Mathlib.Tactic.IntervalCases
import Mathlib.Tactic.Rify
open Finset Filter Topology

-- Important: b should be in ℝ so that it's real div not nat div
lemma geom_sum (n : ℕ) (b : ℝ) (hn : 2 ≤ n) (hb : 2 ≤ b) : ∑ a ∈ Ico 2 n, 1 / (b ^ a) = 1 / (b - 1) - 1 / b - 1 / (b ^ (n - 1) * (b - 1)) := by
  have : ∑ a ∈ Ico 2 n, 1 / (b ^ a) = ∑ a ∈ Ico 2 n, (1 / b) ^ a := by simp
  rw [this, geom_sum_Ico (by grind) hn]
  field_simp
  have : b ^ n ≠ 0 := by positivity
  grind [inv_pow, mul_pow_sub_one]

lemma telescope_sum (n : ℕ) (h : 2 ≤ n) : ∑ a ∈ Ico 2 n, ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ a) = 1 - (1 : ℝ) / (n - 1) - ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) := by
  calc
    ∑ a ∈ Ico 2 n, ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ a) = ∑ b ∈ Ico 2 n, ∑ a ∈ Ico 2 n, (1 : ℝ) / (b ^ a) := sum_comm
    _ = ∑ b ∈ Ico 2 n, ((1 : ℝ) / (b - 1) - (1 : ℝ) / b - (1 : ℝ) / (b ^ (n - 1) * (b - 1))) := by
      apply sum_congr rfl
      exact fun b bico ↦ geom_sum n b h <| Nat.ofNat_le_cast.mpr <| List.left_le_of_mem_range' bico
    _ = ∑ b ∈ Ico 2 n, ((1 : ℝ) / (b - 1) - (1 : ℝ) / b) - ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) := by apply sum_sub_distrib
    _ = 1 - (1 : ℝ) / (n - 1) - ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) := by
      suffices ∑ b ∈ Ico 2 n, ((1 : ℝ) / (b - 1) - (1 : ℝ) / b) = 1 - (1 : ℝ) / (n - 1) by rw [this]
      induction h with
      | refl => norm_num
      | step h ih =>
        rw [sum_Ico_succ_top h, ih]
        simp

-- Exponentials grow quickly
lemma exp_larger (n : ℕ) (h : 2 ≤ n) : (n - 2) * 2 ≤ 2 ^ (n - 1) := by
  induction h <;> grind [mul_pow_sub_one]

-- Another similar lemma
lemma exp_larger' (n : ℕ) (h : 2 ≤ n) : (n - 2) * (n * (n - 1)) < (2 * n - 3) * 2 ^ (n - 1) := by
  by_cases h : 4 < n
  · have : n * (n - 1) * (n - 2) < 2 ^ (n + 1) := by
      induction h
      · decide
      · rename_i m _ _
        suffices (m + 1) * (m * (m - 1)) ≤ 2 * (m - 2) * (m * (m - 1)) by grind
        apply Nat.mul_le_mul_right (m * (m - 1))
        grind
    suffices 2 ^ (n + 1) ≤ (2 * n - 3) * 2 ^ (n - 1) by lia
    have : 2 ^ (n + 1) = 4 * 2 ^ (n - 1) := by grind [mul_pow_sub_one]
    rw [this, mul_le_mul_iff_left₀ (by positivity)]
    grind
  · interval_cases n <;> decide

theorem double_sum : Tendsto (fun n : ℕ ↦ ∑ a ∈ Ico 2 n, ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ a)) atTop (𝓝 (1 : ℝ)) := by
  rw [atTop_basis.tendsto_iff (nhds_basis_Ioo_pos 1)]
  intro ε εpos
  simp only [true_and, Set.mem_Ici, Set.mem_Ioo]
  use max (Nat.ceil (3 / ε)) 2
  intro n nlarge
  have h₁ := le_of_max_le_right nlarge
  rw [telescope_sum n h₁]
  constructor
  · -- Greater than 1 - ε
    suffices (∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) ≤ (n - 2) / 2 ^ (n - 1)) ∧ 1 / (n - 1) + (n - 2) / 2 ^ (n - 1) < ε by grind
    constructor
    · have h₂ b (bico : b ∈ Ico 2 n) : (1 : ℝ) / (b ^ (n - 1) * (b - 1)) ≤ 1 / 2 ^ (n - 1) := by
        rw [← one_div_mul_one_div]
        have h₃ : 2 ≤ (b : ℝ) := by
          norm_cast
          exact (List.left_le_of_mem_range' bico)
        have h₄ : 1 / ((b : ℝ) - 1) ≤ 1 := div_le_one₀ (by lia) |>.mpr (by lia)
        have h₅ : 1 / (b : ℝ) ^ (n - 1) ≤ 1 / 2 ^ (n - 1) := by
          simp_rw [one_div]
          exact inv_anti₀ (by norm_num) <| pow_le_pow_left₀ (by norm_num) h₃ (n - 1)
        grw [h₄, h₅]
        simp
      grw [sum_le_sum h₂]
      simp
      norm_cast
    · by_cases 3 / 2 < ε
      · have : 2 ≤ (n : ℝ) := Nat.ofNat_le_cast.mpr h₁
        have : 1 / ((n : ℝ) - 1) ≤ 1 := div_le_one₀ (by lia) |>.mpr (by lia)
        suffices ((n : ℝ) - 2) / 2 ^ (n - 1) ≤ 1 / 2 by grind
        have := exp_larger n h₁
        rw [← one_mul (2 ^ (n - 1))] at this
        exact div_le_div_iff₀ (by norm_num) (by norm_num) |>.mpr (by norm_cast)
      · have := (div_le_comm₀ (by positivity) εpos).mpr <| Nat.ceil_le.mp (le_of_max_le_left nlarge)
        grw [← this]
        have : 3 / (n : ℝ) - 1 / (n - 1) = (2 * n - 3) / (n * (n - 1)) := by
          field_simp
          grind
        rw [← lt_sub_iff_add_lt', this]
        have h₂ : ((n : ℝ) - 2) * (n * (n - 1)) < (2 * n - 3) * 2 ^ (n - 1) := by
          have := exp_larger' n h₁
          rify at this
          grind [Nat.cast_sub]
        have h₃ : 0 < (2 : ℝ) ^ (n - 1) := by norm_num
        exact div_lt_div_iff₀ h₃ (by nlinarith) |>.mpr h₂
  · -- Less than 1
    have h₂ b (bico : b ∈ Ico 2 n) : 0 ≤ (1 : ℝ) / (b ^ (n - 1) * (b - 1)) := by
      rw [one_div_nonneg]
      have : 1 < (b : ℝ) := Nat.one_lt_cast.mpr (List.left_le_of_mem_range' bico)
      exact mul_nonneg (pow_nonneg (by lia) (n - 1)) (by lia)
    have : 0 < (1 : ℝ) / (n - 1) := by simp [Nat.one_lt_cast.mpr h₁]
    grind [sum_nonneg h₂]