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
102. 102
103. 103
104. 104
105. 105
106. 106
107. 107
108. 108
109. 109
110. 110
111. 111
112. 112
113. 113
114. 114
115. 115
116. 116
117. 117
118. 118
119. 119
120. 120
121. 121
122. 122
123. 123
124. 124
125. 125
126. 126
127. 127
128. 128
129. 129
130. 130

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
  simp
  have : b ^ n ≠ 0 := by
    have : 0 < b ^ n := pow_pos (by linarith) n
    linarith
  grind [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 Finset.sum_congr rfl
      intro b bico
      exact 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 =>
        simp
        norm_num
      | step h ih =>
        rw [sum_Ico_succ_top, ih]
        simp
        exact h

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

-- 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
      · trivial
      · 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 linarith
    have : 2 ^ (n + 1) = 4 * 2 ^ (n - 1) := by grind [mul_pow_sub_one]
    simp [this]
    omega
  · interval_cases n <;> simp

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
  use max (Nat.ceil (3 / ε)) 2
  intro n nlarge
  have h₁ : 2 ≤ n := le_of_max_le_right nlarge
  constructor <;> field_simp <;> rw [telescope_sum n h₁] <;> ring_nf
  · -- Greater than 1 - ε
    have : ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) ≤ (n - 2) / 2 ^ (n - 1) := by
      have h₂ b (bico : b ∈ Ico 2 n) : (1 : ℝ) / (b ^ (n - 1) * (b - 1)) ≤ 1 / 2 ^ (n - 1) := by
        have h₃ : 2 ≤ (b : ℝ) := by
          norm_num
          exact (List.left_le_of_mem_range' bico)
        have h₄ : 1 / ((b : ℝ) - 1) ≤ 1 := (div_le_one₀ (by linarith)).mpr (by linarith)
        have h₅ : 1 / (b : ℝ) ^ (n - 1) ≤ 1 / 2 ^ (n - 1) := by
          simp only [← one_div_pow]
          apply pow_le_pow_left₀ (by norm_num)
          simp
          exact inv_anti₀ (by norm_num) h₃
        rw [← one_div_mul_one_div, mul_comm]
        nth_rw 5 [← one_mul 1]
        rw [mul_div_assoc]
        exact mul_le_mul h₄ h₅ (by norm_num) (by norm_num)
      have := sum_le_sum h₂
      simp
      simp [Nat.cast_sub h₁] at this
      exact this
    ring_nf at this
    have : 1 < 1 + ε - 1 / (n - 1) - (n - 2) / 2 ^ (n - 1) := by
      suffices 1 / (n - 1) + (n - 2) / 2 ^ (n - 1) < ε by linarith
      by_cases hε : 3 / 2 < ε
      · have : 2 ≤ (n : ℝ) := Nat.ofNat_le_cast.mpr h₁
        have : 1 / ((n : ℝ) - 1) ≤ 1 := div_le_one₀ (by linarith) |>.mpr (by linarith)
        suffices ((n : ℝ) - 2) / 2 ^ (n - 1) ≤ 1 / 2 by linarith
        have := exp_larger n h₁
        rw [mul_comm, ← one_mul (2 ^ (n - 1))] at this
        exact (div_le_div_iff₀ (by norm_num) (by norm_num)).mpr (by norm_cast)
      · simp at hε
        have : 3 / n ≤ ε := by
          have : 0 < (n : ℝ) := by
            norm_num
            omega
          exact (div_le_comm₀ this εpos).mpr <| Nat.ceil_le.mp (le_of_max_le_left nlarge)
        suffices 1 / ((n : ℝ) - 1) + ((n : ℝ) - 2) / 2 ^ (n - 1) < 3 / n by linarith
        rw [← lt_sub_iff_add_lt']
        have : 3 / (n : ℝ) - 1 / (n - 1) = (2 * n - 3) / (n * (n - 1)) := by
          have : (n : ℝ) - 1 ≠ 0 := by grind
          field_simp
          ring
        rw [this]
        have h₂ : ((n : ℝ) - 2) * (n * (n - 1)) < (2 * n - 3) * 2 ^ (n - 1) := by
          have h₃ : ↑(2 * n - 3) = 2 * (n : ℝ) - 3 := by
            rw [Nat.cast_sub (by linarith)]
            simp
          have := exp_larger' n h₁
          rify [h₁, h₃] 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₂
    ring_nf at this
    linarith
  · -- Less than 1
    have h₂ b (bico : b ∈ Ico 2 n) : 0 < (1 : ℝ) / (b ^ (n - 1) * (b - 1)) := by
      rw [one_div_pos]
      have : 1 < (b : ℝ) := Nat.one_lt_cast.mpr (List.left_le_of_mem_range' bico)
      exact mul_pos (pow_pos (by linarith) (n - 1)) (by linarith)
    have := sum_nonneg <| fun b a => le_of_lt (h₂ b a)
    have : 0 < (1 : ℝ) / (n - 1) := by simp [Nat.one_lt_cast.mpr h₁]
    grind