/Sum.lean at e889a24a243dcf1fd6e1652edac882f9bfc70b6c

import Mathlib
open Finset BigOperators

-- Pretty standard epsilon delta thing
def converges_to (s : ℕ → ℝ) (a : ℝ) := ∀ ε > 0, ∃ N, ∀ n ≥ N, |s n - a| < ε

-- 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 : converges_to (fun N : ℕ => ∑ a ∈ Ico 2 N, ∑ b ∈ Ico 2 N, 1 / (b ^ a)) 1 := by
  intro ε εpos
  use max (Nat.ceil (3 / ε)) 2
  intro n nlarge
  simp
  have h₁ : 2 ≤ n := le_of_max_le_right nlarge
  rw [abs_lt]
  constructor
  -- Greater than 1 - ε
  field_simp
  rw [telescope_sum n h₁]
  ring_nf
  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₄
  ring_nf at 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 := by
        refine (div_le_one₀ ?_).mpr (by linarith)
        linarith
      suffices ((n : ℝ) - 2) / 2 ^ (n - 1) ≤ 1 / 2 by linarith
      have h₄ : 2 * (n - 2) ≤ 2 ^ (n - 1) := by
        if h₅ : n = 2 then
          simp [h₅]
        else
          exact exp_larger n (by omega)
      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
        if h₅ : n = 2 then
          simp [h₅]
        else if h₆ : n = 3 then
          simp [h₆]
        else if h₇ : n = 4 then
          simp [h₇]
        else
          exact exp_larger' n (by omega)
      apply lt_sub_iff_add_lt'.mp
      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 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