/Sum.lean at ebe50f5e051a8efd668baa2f62241edcd4f9aba6

import Mathlib
open Finset BigOperators Filter Topology

-- Important: b should be in ℝ so that it's real div not nat div
lemma geom_sum (b : ℝ) (N : ℕ) (h₁ : 2 ≤ N) (h₂ : 2 ≤ b) : ∑ a ∈ Ico 2 N, 1 / (b ^ a) = 1 / (b - 1) - 1 / b - 1 / (b ^ (N - 1) * (b - 1)) :=
  calc
    ∑ a ∈ Ico 2 N, 1 / (b ^ a) = ∑ a ∈ Ico 2 N, (1 / b) ^ a := by simp
    _ = ((1 / b) ^ 2 - (1 / b) ^ N) / (1 - (1 / b)) := by
      rw [geom_sum_Ico']
      grind
      exact h₁
    _ = 1 / (b - 1) - 1 / b - b / (b ^ N * (b - 1)) := by
      have : b - 1 ≠ 0 := by linarith
      field_simp
      ring
    _ = 1 / (b - 1) - 1 / b - 1 / (b ^ (N - 1) * (b - 1)) := by
      have h₃ : b ≠ 0 := by linarith
      have h₄ : b ^ N * (b - 1) = b * (b ^ (N - 1) * (b - 1)) := by
        rw [← mul_assoc, ← pow_succ' b (N - 1), Nat.sub_one_add_one]
        linarith
      simp [h₄, div_mul_cancel_left₀ h₃]

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
      have h₁ : 2 ≤ (b : ℝ) := Nat.ofNat_le_cast.mpr (List.left_le_of_mem_range' bico)
      exact geom_sum b N h h₁
    _ = ∑ 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 =>
        field_simp
        linarith
      | step h ih =>
        rw [sum_Ico_succ_top, ih]
        field_simp
        exact h

-- Exponentials grow quickly
lemma exp_larger (n : ℕ) (h : 3 ≤ n) : 2 * (n - 2) ≤ 2 ^ (n - 1) := by
  induction h with
  | refl =>
    trivial
  | step a a_ih =>
    rename_i m
    simp
    simp at a
    rw [← mul_pow_sub_one]
    grind
    omega

-- Another similar lemma
lemma exp_larger' (n : ℕ) (h : 5 ≤ n) : n * (n - 1) * (n - 2) < (2 * n - 3) * 2 ^ (n - 1) := by
  have h₁ : n * (n - 1) * (n - 2) < 2 ^ (n + 1) := by
    induction h with
    | refl =>
      trivial
    | step a a_ih =>
      rename_i m
      simp
      simp at a
      suffices (m + 1) * m * (m - 1) ≤ 2 * (m * (m - 1) * (m - 2)) by grind
      rw [← mul_assoc, mul_assoc]
      nth_rw 3 [mul_comm]
      nth_rw 2 [← mul_assoc]
      apply Nat.mul_le_mul_right (m * (m - 1))
      omega
  suffices 2 ^ (n + 1) ≤ (2 * n - 3) * 2 ^ (n - 1) by linarith
  have h₂ : 2 ^ (n + 1) = 2 * (2 * 2 ^ (n - 1)) := by grind [mul_pow_sub_one]
  rw [h₂, ← mul_assoc]
  apply Nat.mul_le_mul_right (2 ^ (n - 1))
  omega

theorem double_sum : Tendsto (fun N : ℕ => ∑ a ∈ Ico 2 N, ∑ b ∈ Ico 2 N, (1 : ℝ) / (b ^ a)) atTop (𝓝 (1 : ℝ)) := by
  have : atTop.HasBasis (fun _ : ℕ ↦ True) Set.Ici := atTop_basis
  rw [this.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
  · -- Greater than 1 - ε
    field_simp
    rw [telescope_sum n h₁]
    have h₂ : ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) ≤ (n - 2) / 2 ^ (n - 1) := by
      have h₃ : ∀ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) ≤ 1 / 2 ^ (n - 1) := by
        intro b bico
        have h₄ : 2 ≤ (b : ℝ) := by
          norm_num
          exact (List.left_le_of_mem_range' bico)
        have h₅ : 1 / ((b : ℝ) - 1) ≤ 1 := by
          refine (div_le_one₀ ?_).mpr (by linarith)
          linarith
        have h₆ : 1 / (b : ℝ) ^ (n - 1) ≤ 1 / 2 ^ (n - 1) := by
          rw [← one_div_pow, ← 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 h₄ : ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) ≤ ∑ b ∈ Ico 2 n, 1 / 2 ^ (n - 1) := sum_le_sum h₃
      field_simp at h₄
      field_simp
      exact h₄
    have h₃ : 1 < 1 + ε - 1 / (n - 1) - (n - 2) / 2 ^ (n - 1) := by
      suffices 1 / (n - 1) + (n - 2) / 2 ^ (n - 1) < ε by linarith
      if hε : 3 / 2 < ε then
        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 h₄ : 2 * (n - 2) ≤ 2 ^ (n - 1) := by
          by_cases h₅ : n > 2
          · exact exp_larger n h₅
          · grind
        rw [mul_comm, ← one_mul (2 ^ (n - 1))] at h₄
        exact (div_le_div_iff₀ (by norm_num) (by norm_num)).mpr (by norm_cast)
      else
        simp at hε
        have h₃ : 3 / n ≤ ε := by
          have h₄ : 3 / ε ≤ n := Nat.ceil_le.mp (le_of_max_le_left nlarge)
          refine (div_le_comm₀ ?_ εpos).mpr h₄
          norm_num
          omega
        suffices 1 / ((n : ℝ) - 1) + ((n : ℝ) - 2) / 2 ^ (n - 1) < 3 / n by linarith
        have h₄ : n * (n - 1) * (n - 2) < (2 * n - 3) * 2 ^ (n - 1) := by
          by_cases h₅ : n > 4
          · exact exp_larger' n h₅
          · interval_cases n <;> simp
        rw [← lt_sub_iff_add_lt']
        have h₅ : 3 / (n : ℝ) - 1 / (n - 1) = (2 * n - 3) / (n * (n - 1)) := by
          have : (n : ℝ) ≠ 0 := by grind [Nat.cast_ne_zero]
          have : (n : ℝ) - 1 ≠ 0 := by grind
          field_simp
          ring
        rw [h₅]
        have h₆ : ((n : ℝ) - 2) * (n * (n - 1)) < (2 * n - 3) * 2 ^ (n - 1) := by
          nth_rw 1 [mul_comm] at h₄
          rify at h₄
          have h₇ : ↑(n - 1) = (n : ℝ) - 1 := by
            rw [Nat.cast_sub (by linarith)]
            simp
          have h₈ : ↑(2 * n - 3) = 2 * (n : ℝ) - 3 := by
            rw [Nat.cast_sub (by linarith)]
            simp
          have h₉ : ↑(n - 2) = (n : ℝ) - 2 := by norm_cast
          rwa [h₇, h₈, h₉] at h₄
        have h₇ : 0 < (2 : ℝ) ^ (n - 1) := by norm_num
        exact (div_lt_div_iff₀ h₇ (by nlinarith)).mpr h₆
    ring_nf
    ring_nf at h₂
    ring_nf at h₃
    linarith
  · -- Less than 1
    field_simp
    rw [telescope_sum n h₁]
    ring_nf
    field_simp
    have : 0 < (1 : ℝ) / (-1 + n) := by
      field_simp
      omega
    have h₂ : ∀ b ∈ Ico 2 n, 0 < (1 : ℝ) / (b ^ (n - 1) * (b - 1)) := by
      intro b bico
      field_simp
      have : 1 < (b : ℝ) := Nat.one_lt_cast.mpr (List.left_le_of_mem_range' bico)
      have h₃ : 0 < (b : ℝ) ^ (n - 1) := by apply pow_pos (by linarith)
      exact mul_pos h₃ (by linarith)
    have h₃ : 0 ≤ ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) := sum_nonneg (fun b a => le_of_lt (h₂ b a))
    field_simp at h₃
    grind