/Sum.lean at 4bbd02ce5c3a26cdad3fa0b2c354a7304ba259a0

import Mathlib
open Finset BigOperators

-- OK let's do some warm-up
example (b : ℝ) (h : b ≥ 2) : 1 / ((b - 1) * b) = 1 / (b - 1) - 1 / b := by
  have h₁ : b - 1 ≠ 0 := by
    linarith
  field_simp

-- More easy stuff
example (n : ℕ) : (∑ k ∈ Finset.range n, (2 * k + 1)) = n ^ 2 := by
  induction n with
  | zero => trivial
  | succ n ih =>
    rw [sum_range_succ, ih]
    ring

-- How to rw the inside of a sum
example (n : ℕ) : (∑ k ∈ Finset.range n, (2 * k + 1)) = (∑ k ∈ Finset.range n, (2 * (k + 1) - 1)) := by
  apply Finset.sum_congr rfl
  intro k hk
  omega

-- 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
-- This was so excruciating wow
lemma geom_sum (b : ℝ) (N : ℕ) (h₁ : N ≥ 2) (h₂ : b ≥ 2) : ∑ 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']
      by_contra h
      simp at h
      linarith
      exact h₁
    _ = 1 / (b - 1) - 1 / b - b / (b ^ N * (b - 1)) := by
      have h₃ : 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 : N ≥ 2): ∑ 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₁ : (b : ℝ) ≥ 2 := 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)) := sum_sub_distrib
    _ = 1 - (1 : ℝ) / (N - 1) - ∑ b ∈ Ico 2 N, (1 : ℝ) / (b ^ (N - 1) * (b - 1)) := by
      have h₁ : ∑ b ∈ Ico 2 N, ((1 : ℝ) / (b - 1) - (1 : ℝ) / b) = 1 - (1 : ℝ) / (N - 1) := by
        induction h with
        | refl =>
          field_simp
          linarith
        | step h ih =>
          rw [sum_Ico_succ_top, ih]
          field_simp
          exact h
      rw [h₁]

-- Exponentials grow quickly
lemma exp_larger (n : ℕ) (h : n ≥ 3) : 2 * (n - 2) ≤ 2 ^ (n - 1) := by
  induction h with
  | refl =>
    trivial
  | step a a_ih =>
    rename_i m
    simp
    simp at a
    have h₁ : 2 * 2 ^ (m - 1) = 2 ^ m := by
      apply mul_pow_sub_one
      linarith
    rw [←h₁]
    apply Nat.mul_le_mul_left 2
    have h₂ : m - 1 ≤ 2 * (m - 2) := by
      omega
    linarith

-- Another similar lemma
lemma exp_larger' (n : ℕ) (h : n ≥ 5) : 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
      have h₂ : (m + 1) * m * (m - 1) ≤ 2 * (m * (m - 1) * (m - 2)) := by
        rw [←mul_assoc]
        rw [mul_assoc]
        nth_rw 3 [mul_comm]
        nth_rw 2 [←mul_assoc]
        apply Nat.mul_le_mul_right (m * (m - 1))
        omega
      have h₃ : 2 ^ (m + 2) = 2 * 2 ^ (m + 1) := Nat.pow_succ'
      linarith
  have h₂ : 2 ^ (n + 1) ≤ (2 * n - 3) * 2 ^ (n - 1) := by
    have h₃ : 2 * (2 * 2 ^ (n - 1)) = 2 ^ (n + 1) := by
      rw [mul_pow_sub_one]
      omega
      omega
    rw [←h₃, ←mul_assoc]
    apply Nat.mul_le_mul_right (2 ^ (n - 1))
    omega
  linarith

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₁ : n ≥ 2 := by
    exact 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₄ : b ≥ 2 := List.left_le_of_mem_range' bico
      have h₅ : (b : ℝ) ≥ 2 := by
        norm_num
        linarith
      have h₆ : 1 / ((b : ℝ) - 1) ≤ 1 := by
        have h₇ : (b : ℝ) - 1 ≥ 1 := by
          linarith
        refine (div_le_one₀ ?_).mpr h₇
        linarith
      have h₇ : 1 / (b : ℝ) ^ (n - 1) ≤ 1 / 2 ^ (n - 1) := by
        rw [←one_div_pow, ←one_div_pow]
        apply pow_le_pow_left₀
        norm_num
        simp
        apply inv_anti₀ at h₅
        exact h₅
        norm_num
      have h₈ : 0 ≤ 1 / (b : ℝ) ^ (n - 1) := by
        norm_num
      have h₉ : (0 : ℝ) ≤ 1 := by
        norm_num
      rw [←one_div_mul_one_div, mul_comm]
      nth_rw 5 [←one_mul 1]
      rw [mul_div_assoc]
      exact mul_le_mul h₆ h₇ h₈ h₉
    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
    have h₄ : 1 / (n - 1) + (n - 2) / 2 ^ (n - 1) < ε := by
      if hε : ε > 3 / 2 then
        have h₅ : (n : ℝ) ≥ 2 := by
          norm_num
          linarith
        have h₆ : 1 / ((n : ℝ) - 1) ≤ 1 := by
          have h₇ : (n : ℝ) - 1 ≥ 1 := by
            linarith
          refine (div_le_one₀ ?_).mpr h₇
          linarith
        have h₇ : ((n : ℝ) - 2) / 2 ^ (n - 1) ≤ 1 / 2 := by
          have h₈ : 2 * (n - 2) ≤ 2 ^ (n - 1) := by
            if h₉ : n = 2 then
              rw [h₉]
              norm_num
            else
              have h₁₀ : 3 ≤ n := by
                omega
              exact exp_larger n h₁₀
          rw [mul_comm, ←one_mul (2 ^ (n - 1))] at h₈
          have h₉ : 0 < (2 : ℝ) ^ (n - 1) := by
            norm_num
          have h₁₀ : 0 < (2 : ℝ) := by
            norm_num
          have h₁₁ : ((n : ℝ) - 2) * 2 ≤ 1 * 2 ^ (n - 1) := by
            norm_cast
          exact (div_le_div_iff₀ h₉ h₁₀).mpr h₁₁
        linarith
      else
        simp at hε
        have h₃ : n ≥ Nat.ceil (3 / ε) := by
          exact le_of_max_le_left nlarge
        have h₄ : 3 / n ≤ ε := by
          have h₅ : 3 / ε ≤ n := by
            exact Nat.ceil_le.mp h₃
          refine (div_le_comm₀ ?_ εpos).mpr h₅
          norm_num
          linarith
        have h₅ : 1 / ((n : ℝ) - 1) + ((n : ℝ) - 2) / 2 ^ (n - 1) < 3 / n := by
          have h₆ : n * (n - 1) * (n - 2) < (2 * n - 3) * 2 ^ (n - 1) := by
            if h₇ : n = 2 then
              rw [h₇]
              norm_num
            else if h₈ : n = 3 then
              rw [h₈]
              norm_num
            else if h₉ : n = 4 then
              rw [h₉]
              norm_num
            else
              have h₁₀ : 5 ≤ n := by
                omega
              exact exp_larger' n h₁₀
          apply lt_sub_iff_add_lt'.mp
          have h₇ : 3 / (n : ℝ) - 1 / (n - 1) = (2 * n - 3) / (n * (n - 1)) := by
            have h₈ : (n : ℝ) ≠ 0 := by
              rw [Nat.cast_ne_zero]
              linarith
            have h₉ : (n : ℝ) - 1 ≠ 0 := by
              have h₁₀ : 0 < n := by
                linarith
              have h₁₁ : ↑(n - 1) = (n : ℝ) - 1 := by
                rw [Nat.cast_sub h₁₀]
                simp
              rw [←h₁₁]
              norm_cast
              omega
            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
              have h₁₀ : 0 < n := by
                linarith
              rw [Nat.cast_sub h₁₀]
              simp
            have h₁₀ : ↑(2 * n - 3) = 2 * (n : ℝ) - 3 := by
              have h₁₁ : 2 < 2 * n := by
                linarith
              rw [Nat.cast_sub h₁₁]
              simp
            have h₁₁ : ↑(n - 2) = (n : ℝ) - 2 := by
              norm_cast
            rw [h₉, h₁₀, h₁₁] at h₆
            exact h₆
          have h₉ : 0 < (2 : ℝ) ^ (n - 1) := by
            norm_num
          have h₁₀ :  0 < (n : ℝ) * (n - 1) := by
            nlinarith
          exact (div_lt_div_iff₀ h₉ h₁₀).mpr h₈
        linarith
    linarith
  ring_nf at h₃
  linarith
  -- Less than 1
  field_simp
  rw [telescope_sum n h₁]
  ring_nf
  field_simp
  have h₂ : (1 : ℝ) / (-1 + n) > 0 := by
    field_simp
    linarith
  have h₃ : ∀ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) > 0 := by
    intro b bico
    field_simp
    have h₄ : b ≥ 2 := List.left_le_of_mem_range' bico
    have h₅ : (b : ℝ) > 1 := Nat.one_lt_cast.mpr h₄
    have h₆ : (b : ℝ) ^ (n - 1) > 0 := by
      have h₇ : (b : ℝ) > 0 := by
        linarith
      apply pow_pos h₇
    have h₇ : (b : ℝ) - 1 > 0 := by
      linarith
    apply mul_pos h₆ h₇
  have h₄ : ∑ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) ≥ 0 := by
    have h₅ : ∀ b ∈ Ico 2 n, (1 : ℝ) / (b ^ (n - 1) * (b - 1)) ≥ 0 := fun b a => le_of_lt (h₃ b a)
    exact sum_nonneg h₅
  ring_nf at h₄
  field_simp at h₄
  have h₅ : -(1 : ℝ) / (-1 + n) < 0 := by
    nth_rw 1 [← mul_neg_one, mul_comm, mul_div_assoc]
    linarith
  linarith