Erdos481.lean

import Mathlib.Data.Nat.Basic
import Mathlib.Tactic
import Mathlib.Data.Real.Basic
import Mathlib.Algebra.Order.Field.Basic
import Mathlib.Algebra.Order.Field.Rat
import Mathlib.Analysis.SpecialFunctions.Exp
import Mathlib.Analysis.SpecialFunctions.Pow.Asymptotics
import Mathlib.Analysis.SpecialFunctions.PolynomialExp
import Mathlib.Analysis.SpecificLimits.Normed
import Mathlib.Order.Filter.AtTopBot.Basic
import Mathlib.Analysis.SpecialFunctions.Pow.Real
import Mathlib.NumberTheory.Harmonic.Bounds
import Mathlib.Data.Nat.Cast.Field
import Mathlib.Analysis.SpecialFunctions.Log.Basic
import Mathlib.Data.Rat.Defs
import Mathlib.Algebra.BigOperators.Group.List.Basic
import Mathlib.Algebra.Group.Defs
import Mathlib.Data.List.Basic
import Mathlib.Control.Basic
import Mathlib.Data.List.FinRange
import Mathlib.Algebra.BigOperators.Group.Finset.Defs
import Mathlib.Data.Fintype.Basic
import Mathlib.Data.Fin.Basic
import Mathlib.Algebra.BigOperators.Fin
import Mathlib.Algebra.BigOperators.Ring.List
import Mathlib.Data.List.Range
import Mathlib.Algebra.BigOperators.Group.Finset.Basic
import Mathlib.Data.Nat.Cast.Basic
import Mathlib.Algebra.Order.Ring.Nat
import Mathlib.Analysis.SpecificLimits.Basic
import Mathlib.Data.List.Sort
import Mathlib.Data.List.NodupEquivFin
import Mathlib.Algebra.Order.BigOperators.Group.Finset
import Mathlib.Data.Finset.Sort
import Mathlib.Algebra.Field.Rat
import Mathlib.Algebra.Order.Ring.Rat
import Mathlib.Data.List.Permutation
import Mathlib.Data.List.OfFn
import Mathlib.Data.Rat.Cast.Order
import Mathlib.Algebra.Order.BigOperators.Group.List
import Mathlib.Algebra.Order.Archimedean.Basic
import Mathlib.Data.Rat.Cast.Defs
import Mathlib.Data.Rat.BigOperators
import Mathlib.Algebra.Field.Defs
import Mathlib.Algebra.Ring.Rat
import Mathlib.Data.Finset.Basic
import Mathlib.Analysis.PSeries
import Mathlib.Data.List.Flatten

set_option linter.style.longLine false
set_option linter.unusedVariables false
set_option linter.unnecessarySimpa false
set_option linter.unusedSimpArgs false

open BigOperators

theorem child_ge_of_parent_ge (a b x M : ℕ) (ha : a ≥ 1) (hx : x ≥ M) : a * x + b ≥ M := by
  nlinarith

theorem div_sub_div_sq_ge_iff_mul (σ ρ C n : ℝ) (hn : 0 < n) :
  σ / n - C / n ^ 2 ≥ ρ / n ↔ n * (σ - ρ) ≥ C := by
  have hn0 : (n : ℝ) ≠ 0 := ne_of_gt hn
  -- useful algebraic identity: σ*n - ρ*n = n*(σ - ρ)
  have hσρ : σ * n - ρ * n = n * (σ - ρ) := by
    ring
  constructor
  · intro h
    -- rewrite to ≤ form
    have h1 : ρ / n ≤ σ / n - C / n ^ 2 := by
      simpa [ge_iff_le] using h
    -- clear denominators to get polynomial inequality
    have h2 : ρ * n ≤ σ * n - C := by
      have h1' := h1
      field_simp [hn0] at h1'
      simpa [pow_two, mul_comm, mul_left_comm, mul_assoc, sub_eq_add_neg] using h1'
    -- From ρ*n ≤ σ*n - C, get ρ*n + C ≤ σ*n
    have hA : ρ * n + C ≤ σ * n := (le_sub_iff_add_le).1 h2
    -- Rearrange sum
    have hA' : C + ρ * n ≤ σ * n := by
      simpa [add_comm, add_left_comm, add_assoc] using hA
    -- From C + ρ*n ≤ σ*n, get C ≤ σ*n - ρ*n
    have h3 : C ≤ σ * n - ρ * n := (le_sub_iff_add_le).2 hA'
    -- Replace σ*n - ρ*n with n*(σ - ρ)
    have h4 : C ≤ n * (σ - ρ) := by
      simpa [hσρ] using h3
    -- Back to ≥ form
    simpa [ge_iff_le] using h4
  · intro h
    -- rewrite to ≤ form
    have h1 : C ≤ n * (σ - ρ) := by
      simpa [ge_iff_le] using h
    -- rewrite RHS using the algebraic identity
    have h1a : C ≤ σ * n - ρ * n := by
      simpa [hσρ] using h1
    -- from C ≤ σ*n - ρ*n, get C + ρ*n ≤ σ*n
    have hA : C + ρ * n ≤ σ * n := (le_sub_iff_add_le).1 h1a
    -- rearrange sum
    have hA' : ρ * n + C ≤ σ * n := by
      simpa [add_comm, add_left_comm, add_assoc] using hA
    -- from ρ*n + C ≤ σ*n, get ρ*n ≤ σ*n - C
    have hpoly : ρ * n ≤ σ * n - C := (le_sub_iff_add_le).2 hA'
    -- turn polynomial inequality into the fractional one
    have h3 : ρ / n ≤ σ / n - C / n ^ 2 := by
      field_simp [hn0]
      -- goal is now ρ * n ≤ σ * n - C
      simpa [pow_two, mul_comm, mul_left_comm, mul_assoc, sub_eq_add_neg] using hpoly
    -- back to ≥ form
    simpa [ge_iff_le] using h3

theorem exp_dominates_linear (C : ℝ) (hC : C > 0) (ρ : ℝ) (hρ : ρ > 1) (A B : ℝ) : ∃ k : ℕ, C * ρ ^ k > A * (k : ℝ) + B := by
  have h_rho_pos : 0 < ρ := zero_lt_one.trans hρ

  -- Limit 1: k / ρ^k -> 0
  have lim1 : Filter.Tendsto (fun k : ℕ => (k : ℝ) / ρ ^ k) Filter.atTop (nhds 0) := by
    convert tendsto_pow_const_div_const_pow_of_one_lt 1 hρ using 1
    ext k
    simp

  -- Limit 2: (ρ⁻¹)^k -> 0
  have lim2 : Filter.Tendsto (fun k : ℕ => (ρ⁻¹) ^ k) Filter.atTop (nhds 0) := by
    apply tendsto_pow_atTop_nhds_zero_of_abs_lt_one
    rw [abs_inv, abs_of_pos h_rho_pos]
    exact inv_lt_one_of_one_lt₀ hρ

  -- Define f
  let f := fun (k : ℕ) => (A * k + B) / (C * ρ ^ k)

  -- Show f -> 0
  have lim_f : Filter.Tendsto f Filter.atTop (nhds 0) := by
    let g := fun (k : ℕ) => (A / C) * ((k : ℝ) / ρ ^ k) + (B / C) * (ρ⁻¹) ^ k
    have h_eq : ∀ k, f k = g k := by
      intro k
      dsimp [f, g]
      rw [inv_pow]
      field_simp [pow_ne_zero k (ne_of_gt h_rho_pos), ne_of_gt hC]
    rw [show (0 : ℝ) = (A / C) * 0 + (B / C) * 0 by ring]
    apply Filter.Tendsto.congr (f₁ := g)
    · intro k; rw [h_eq k]
    · apply Filter.Tendsto.add
      · apply Filter.Tendsto.const_mul; exact lim1
      · apply Filter.Tendsto.const_mul; exact lim2

  -- Extract k
  have : Filter.Eventually (fun k => f k < 1) Filter.atTop := lim_f (Iio_mem_nhds zero_lt_one)

  rcases this.exists with ⟨k, hk⟩
  use k
  dsimp [f] at hk
  rw [div_lt_one (mul_pos hC (pow_pos h_rho_pos k))] at hk
  exact hk

theorem flatMap_pure_eq_map {α β : Type} (l : List α) (f : α → β) : l.flatMap (fun x => [f x]) = l.map f := by
  induction l with
  | nil =>
    simp [List.flatMap]
  | cons x xs ih =>
    simp [List.flatMap] at ih
    simpa [List.flatMap] using ih

theorem harmonic_Q_le_one_add_log (n : ℕ) (hn : 1 ≤ n) : (harmonic n : ℝ) ≤ 1 + Real.log (n : ℝ) := by
  -- We can directly use the existing bound from Mathlib.
  simpa using harmonic_le_one_add_log n

theorem inv_add_ge_sub_square (a t : ℚ) (ha : 0 < a) (ht : 0 ≤ t) : (1 : ℚ) / (a + t) ≥ 1 / a - t / (a ^ 2) := by
  have ha2 : 0 < a^2 := pow_pos ha 2
  have hat : 0 < a + t := add_pos_of_pos_of_nonneg ha ht
  have hden : 0 < a^2 * (a + t) := mul_pos ha2 hat
  rw [ge_iff_le, ← mul_le_mul_iff_left₀ hden]
  field_simp [ne_of_gt ha, ne_of_gt hat]
  ring_nf
  have : - (t ^ 2) ≤ 0 := neg_nonpos.mpr (sq_nonneg t)
  linarith

theorem inv_add_ge_sub_square_helper (a x b : ℚ) (hx : x > 0) (ha : a > 0) (hb : b ≥ 0) : 1 / (a * x + b) ≥ 1 / (a * x) - b / (a ^ 2 * x ^ 2) := by
  have hax : a * x > 0 := mul_pos ha hx
  have d1 : a * x + b > 0 := add_pos_of_pos_of_nonneg hax hb
  have d2 : a * x ≠ 0 := ne_of_gt hax
  have d3 : a^2 * x^2 > 0 := by nlinarith

  rw [ge_iff_le]

  have key : 1 / (a * x) - b / (a ^ 2 * x ^ 2) = (a * x - b) / (a ^ 2 * x ^ 2) := by
    field_simp [d2, ne_of_gt d3]

  rw [key]

  have denom_pos : (a ^ 2 * x ^ 2) * (a * x + b) > 0 := mul_pos d3 d1
  rw [← mul_le_mul_iff_left₀ denom_pos]

  field_simp [ne_of_gt d3, ne_of_gt d1]
  ring_nf
  nlinarith

theorem list_bind_sum_map {α β γ : Type} [AddCommMonoid γ] (L : List α) (f : α → List β) (g : β → γ) :
  ((L >>= f).map g).sum = (L.map (fun a => ((f a).map g).sum)).sum := by
  induction L with
  | nil =>
      simp
  | cons a L ih =>
      -- transform ih to use flatMap instead of >>=
      have ih' := ih
      simp [List.bind_eq_flatMap] at ih'
      -- now simplify main goal using ih'
      simpa [List.bind_eq_flatMap, ih']

theorem flatMap_sum_map {α β γ : Type} [AddCommMonoid γ] (L : List α) (f : α → List β) (g : β → γ) :
  ((List.flatMap f L).map g).sum = (L.map (fun a => ((f a).map g).sum)).sum := by
  -- reduce to list_bind_sum_map, rewriting bind as flatMap
  simpa [List.bind_eq_flatMap] using
    (list_bind_sum_map (L := L) (f := f) (g := g))

theorem list_sum_comm {α β γ : Type} [AddCommMonoid γ] (L : List α) (M : List β) (f : α → β → γ) :
  (L.map (fun a => (M.map (fun b => f a b)).sum)).sum = (M.map (fun b => (L.map (fun a => f a b)).sum)).sum := by
  induction L with
  | nil => simp
  | cons a L ih =>
    simp [ih, List.sum_map_add]

theorem list_sum_finRange_eq_finset_sum {β : Type} [AddCommMonoid β] (n : ℕ) (f : Fin n → β) : ((List.finRange n).map f).sum = ∑ i : Fin n, f i := by
  classical
  -- We use that `List.ofFn f` is equal to `(List.finRange n).map f`,
  -- and the lemma `List.sum_ofFn` relating the list sum to the finset sum over `Fin n`.
  -- First, rewrite the goal using `List.ofFn`.
  have h₁ : ((List.finRange n).map f).sum = (List.ofFn f).sum := by
    simpa [List.ofFn_eq_map (f := f)]
  -- Next, use `List.sum_ofFn` to identify `(List.ofFn f).sum` with the finset sum.
  have h₂ : (List.ofFn f).sum = ∑ i : Fin n, f i := by
    simpa using (List.sum_ofFn (f := f))
  -- Combine both equalities.
  -- From `h₁` and `h₂` we get the desired statement by transitivity.
  simpa [h₁] using h₂

lemma my_sum_add_distrib {β : Type} [AddCommMonoid β] {ι : Type} (s : Finset ι) (f g : ι → β) :
  ∑ i ∈ s, (f i + g i) = ∑ i ∈ s, f i + ∑ i ∈ s, g i :=
  Finset.sum_add_distrib

theorem list_sum_finset_sum_comm {α β : Type} [AddCommMonoid β] (L : List α) (r : ℕ) (f : α → Fin r → β) :
  (L.map (fun a => ∑ i : Fin r, f a i)).sum = ∑ i : Fin r, (L.map (fun a => f a i)).sum := by
  induction L with
  | nil => simp
  | cons a L ih =>
    simp [ih]
    exact (my_sum_add_distrib Finset.univ (f a) (fun i => (L.map (fun a => f a i)).sum)).symm

theorem list_sum_map_const_mul {α : Type} (L : List α) (c : ℝ) (f : α → ℝ) : (L.map (fun x => c * f x)).sum = c * (L.map f).sum := by
  induction L with
  | nil => simp
  | cons a L ih =>
    simp
    rw [ih]
    rw [mul_add]

theorem list_sum_range_eq_finset_sum_range {β : Type} [AddCommMonoid β] (f : ℕ → β) (n : ℕ) :
  ((List.range n).map f).sum = ∑ i ∈ Finset.range n, f i := by
  induction n with
  | zero =>
      simp
  | succ n ih =>
      simpa [List.range_succ, ih, Finset.range_add_one, add_comm, add_left_comm, add_assoc]

theorem harmonic_list_sum_Q (n : ℕ) : ((List.range n).map (fun x => (1 : ℚ) / ((x : ℚ) + 1))).sum = harmonic n := by
  classical
  -- We will prove the statement by working directly with the function
  -- `f i = (i+1)⁻¹ : ℚ` on natural numbers, and then rewriting the integrand
  -- at the end.
  let f : ℕ → ℚ := fun i => ((i + 1 : ℚ))⁻¹
  have main : ((List.range n).map f).sum = harmonic n := by
    -- Induction on `n`.
    induction n with
    | zero =>
        -- `List.range 0` is empty and `harmonic 0 = 0`.
        simp [f, harmonic_zero]
    | succ n ih =>
        -- `List.range (n+1)` is `List.range n` with `n` appended.
        have hrange : List.range (n + 1) = List.range n ++ [n] := by
          simpa using List.range_succ (n := n)
        -- Compute the sum at `n+1` in terms of the sum at `n`.
        have hsum : ((List.range (n + 1)).map f).sum
            = ((List.range n).map f).sum + f n := by
          simp [hrange, List.map_append, List.sum_append, f]
        -- Now connect with the harmonic recurrence.
        calc
          ((List.range (n + 1)).map f).sum
              = ((List.range n).map f).sum + f n := hsum
          _ = harmonic n + f n := by simpa [ih]
          _ = harmonic n + ((n + 1 : ℚ))⁻¹ := by rfl
          _ = harmonic (n + 1) := by
                -- `harmonic_succ n : harmonic (n + 1) = harmonic n + (↑(n + 1))⁻¹`.
                simpa using (harmonic_succ n).symm
  -- Finally, rewrite the integrand of the list sum into the desired form.
  -- For each `x : ℕ`, `(x+1 : ℚ)⁻¹ = 1 / ((x : ℚ) + 1)`.
  have hpt (x : ℕ) : f x = (1 : ℚ) / ((x : ℚ) + 1) := by
    simp [f, one_div, add_comm, add_left_comm, add_assoc]
  -- Use `hpt` to rewrite the mapped list; this works elementwise.
  have hmap : (List.range n).map f =
      (List.range n).map (fun x => (1 : ℚ) / ((x : ℚ) + 1)) := by
    -- We show pointwise equality of the mapped elements via induction on the list.
    -- This avoids dealing with the coercions that show up in the pretty-printer.
    clear main
    -- Induct on `n` (equivalently, on `List.range n`).
    induction n with
    | zero =>
        simp
    | succ n ih =>
        have hrange : List.range (n + 1) = List.range n ++ [n] := by
          simpa using List.range_succ (n := n)
        simp [hrange, hpt, ih]
  -- Now put everything together.
  have := main
  simpa [f, hmap, one_div, add_comm, add_left_comm, add_assoc]

theorem log_nat_pow (r k : ℕ) (hr : 2 ≤ r) : Real.log ((r : ℝ) ^ k) = (k : ℝ) * Real.log (r : ℝ) := by
  -- We can directly use the existing lemma `Real.log_pow`
  -- which states: `Real.log (x ^ n) = n * Real.log x` for `x : ℝ`, `n : ℕ`.
  have := Real.log_pow (r : ℝ) k
  -- `this` has type `Real.log ((r : ℝ) ^ k) = (k : ℕ) * Real.log (r : ℝ)`.
  -- Coerce `k` from `ℕ` to `ℝ` on the right-hand side.
  -- We rewrite using `Nat.cast_mul` and `one_mul` etc., but actually
  -- the RHS already matches `(k : ℝ) * Real.log (r : ℝ)` after coercion.
  -- So we can just use `exact_mod_cast` or rewrite explicitly.
  -- However, `Real.log_pow` already has `k` as a natural, and multiplication
  -- is in `ℝ`, so the type matches exactly, and `this` is our goal.
  simpa using this

theorem geometric_growth_contradicts_harmonic_bound (C ρ : ℝ) (hC : C > 0) (hρ : ρ > 1) (r : ℕ) (hr : r ≥ 2) : ∃ n : ℕ, C * ρ ^ n > 1 + Real.log ((r : ℝ) ^ n) := by
  obtain ⟨n, hn⟩ := exp_dominates_linear C hC ρ hρ (Real.log r) 1
  use n
  rw [log_nat_pow r n hr]
  rw [add_comm 1, mul_comm (n : ℝ)]
  exact hn

theorem map_one_div_coe_list (L : List ℕ) :
  (L.map (fun x => (1 : ℚ) / (x : ℚ))) =
    ((L.map (fun x : ℕ => (x : ℚ))).map (fun q : ℚ => (1 : ℚ) / q)) := by
  -- First show that the `do`-notation list is just the coercion map.
  have hL :
      (do let a ← L; pure (↑a : ℚ)) = List.map (fun x : ℕ => (x : ℚ)) L := by
    -- Rewrite the monadic bind as `flatMap`, then as a plain `map`.
    -- `do let a ← L; pure (↑a : ℚ)` desugars to `L >>= fun a => pure (↑a : ℚ)`.
    -- Step 1: turn bind into flatMap.
    have h1 :
        (do let a ← L; pure (↑a : ℚ)) =
          List.flatMap (fun a : ℕ => pure (↑a : ℚ)) L := by
      -- `bind_eq_flatMap` has type `l >>= f = l.flatMap f`.
      -- The `do`-block is syntactic sugar for `>>=`.
      simpa [List.bind_eq_flatMap]  -- no RHS, goal already has desired form
    -- Step 2: identify this `flatMap` with a plain `map`.
    have h2 :
        List.flatMap (fun a : ℕ => pure (↑a : ℚ)) L =
          List.map (fun x : ℕ => (x : ℚ)) L := by
      -- Use `flatMap_pure_eq_map` with `f := (fun x : ℕ => (x : ℚ))`.
      simpa [Function.comp] using
        (List.flatMap_pure_eq_map (f := fun x : ℕ => (x : ℚ)) (l := L))
    exact h1.trans h2
  -- Now apply `List.map` with `q ↦ 1/q` to both sides.
  have h1div :=
    congrArg (List.map (fun q : ℚ => (1 : ℚ) / q)) hL
  -- This is exactly the goal (up to definitional equality).
  simpa using h1div

theorem one_le_pow_of_two_le (r k : ℕ) (hr : 2 ≤ r) : 1 ≤ r ^ k := by
  have h1 : 1 ≤ r := le_trans (by decide : 1 ≤ 2) hr
  induction k with
  | zero =>
      simp
  | succ k hk =>
      have hk' : 1 ≤ r ^ k * r := by
        have hmul : 1 * 1 ≤ r ^ k * r := Nat.mul_le_mul hk h1
        simpa using hmul
      simpa [Nat.pow_succ] using hk'

theorem recurrence_algebraic_bound (S : ℕ → ℝ) (σ D : ℝ) (hσ : 1 < σ) (h_rec : ∀ k, S (k + 1) ≥ σ * S k - D) : ∀ k, S k ≥ σ ^ k * (S 0 - D / (σ - 1)) + D / (σ - 1) := by
  intro k
  induction k with
  | zero =>
    simp
  | succ k ih =>
    specialize h_rec k
    have hσ_pos : 0 < σ := by linarith
    calc
      S (k + 1) ≥ σ * S k - D := h_rec
      _ ≥ σ * (σ ^ k * (S 0 - D / (σ - 1)) + D / (σ - 1)) - D := by nlinarith
      _ = σ ^ (k + 1) * (S 0 - D / (σ - 1)) + σ * (D / (σ - 1)) - D := by ring
      _ = σ ^ (k + 1) * (S 0 - D / (σ - 1)) + D / (σ - 1) := by
        field_simp [sub_ne_zero.mpr (ne_of_gt hσ)]
        ring

theorem recurrence_with_decay_counterexample : ∃ (u : ℕ → ℝ) (σ : ℝ) (D : ℕ → ℝ), 1 < σ ∧ Filter.Tendsto D Filter.atTop (nhds 0) ∧ (∀ n, u n > 0) ∧ (∀ n, u (n + 1) ≥ σ * u n - D n) ∧ ¬ Filter.Tendsto u Filter.atTop Filter.atTop := by
  let u : ℕ → ℝ := fun n => 1 / (n + 1 : ℝ)
  let σ : ℝ := 2
  let D : ℕ → ℝ := fun n => σ * u n - u (n + 1)
  use u, σ, D
  refine ⟨by norm_num, ?_, ?_, ?_, ?_⟩
  · -- Tendsto D
    have h_u : Filter.Tendsto u Filter.atTop (nhds 0) := tendsto_one_div_add_atTop_nhds_zero_nat
    dsimp [D]
    rw [show (fun n ↦ σ * u n - u (n + 1)) = (fun n ↦ σ * u n) - (fun n ↦ u (n + 1)) by ext; simp]
    rw [show (0 : ℝ) = σ * 0 - 0 by simp]
    apply Filter.Tendsto.sub
    · apply Filter.Tendsto.const_mul
      exact h_u
    · -- Fix composition: `fun n ↦ u (n + 1)` is `u ∘ (fun n ↦ n + 1)`
      have : (fun n ↦ u (n + 1)) = u ∘ (fun n ↦ n + 1) := rfl
      rw [this]
      apply Filter.Tendsto.comp h_u
      exact Filter.tendsto_add_atTop_nat 1
  · -- u n > 0
    intro n
    dsimp [u]
    apply one_div_pos.mpr
    exact Nat.cast_add_one_pos n
  · -- Inequality
    intro n
    dsimp [D]
    linarith
  · -- Not tendsto atTop
    intro h_div
    have h_u : Filter.Tendsto u Filter.atTop (nhds 0) := tendsto_one_div_add_atTop_nhds_zero_nat
    rw [Filter.tendsto_atTop_atTop] at h_div
    obtain ⟨N, hN⟩ := h_div 1
    have h2 : ∀ᶠ n in Filter.atTop, dist (u n) 0 < 0.5 := Metric.tendsto_nhds.mp h_u 0.5 (by norm_num)
    rw [Filter.eventually_atTop] at h2
    obtain ⟨M, hM⟩ := h2
    let K := max N M
    have hK1 : 1 ≤ u K := hN K (le_max_left N M)
    have hK2 : dist (u K) 0 < 0.5 := hM K (le_max_right N M)
    rw [Real.dist_0_eq_abs] at hK2
    have : u K > 0 := by dsimp [u]; apply one_div_pos.mpr; exact Nat.cast_add_one_pos K
    rw [abs_of_pos this] at hK2
    linarith

def seqA (r : ℕ) (a b : Fin r → ℕ) : ℕ → List ℕ
| 0 => [1]
| n + 1 => (List.finRange r) >>= (fun i => (seqA r a b n).map (fun x => a i * x + b i))

lemma nat_two_mul (n : ℕ) : 2 * n = n + n := by
  induction n with
  | zero =>
      simp
  | succ n ih =>
      simp [Nat.mul_succ, ih, Nat.add_comm, Nat.add_left_comm, Nat.add_assoc]

lemma nat_mul_two (n : ℕ) : n * 2 = n + n := by
  simpa [Nat.mul_comm] using nat_two_mul n

theorem min_element_growth (r : ℕ) (a b : Fin r → ℕ) (ha : ∀ i, a i ≥ 1) (h_no_id : ∀ i, ¬ (a i = 1 ∧ b i = 0)) (k : ℕ) : ∀ x ∈ seqA r a b k, x ≥ k + 1 := by
  intro x
  induction k generalizing x with
  | zero =>
      intro hx
      simp [seqA] at hx
      subst hx
      simp
  | succ k ih =>
      intro hx
      -- hx : x ∈ seqA r a b (k + 1)
      simp [seqA] at hx
      rcases hx with ⟨i, y, hy, rfl⟩
      -- ih : ∀ x ∈ seqA r a b k, x ≥ k + 1
      have ha_i : 1 ≤ a i := ha i
      have hy_ge : k + 1 ≤ y := ih y hy
      -- so k + 2 ≤ y + 1
      have hy_ge' : k + 2 ≤ y + 1 := Nat.succ_le_succ hy_ge
      -- also y ≥ 1
      have hy_pos : 1 ≤ y := by
        have one_le_k1 : 1 ≤ k + 1 := Nat.succ_le_succ (Nat.zero_le k)
        exact le_trans one_le_k1 hy_ge
      -- show y + 1 ≤ a i * y + b i
      have step_ge : y + 1 ≤ a i * y + b i := by
        by_cases h_ai1 : a i = 1
        · -- case a i = 1, then b i ≥ 1 from h_no_id
          have hb_ne : b i ≠ 0 := by
            intro hb
            apply h_no_id i
            exact And.intro h_ai1 hb
          have hb_pos : 0 < b i := Nat.pos_of_ne_zero hb_ne
          have hb_ge1 : 1 ≤ b i := Nat.succ_le_of_lt hb_pos
          have h1 : y + 1 ≤ y + b i := by
            exact Nat.add_le_add_left hb_ge1 y
          have : y + 1 ≤ a i * y + b i := by
            simpa [h_ai1] using h1
          exact this
        · -- case a i ≠ 1, then a i ≥ 2
          have hgt : 1 < a i := lt_of_le_of_ne ha_i (Ne.symm h_ai1)
          have ha2 : 2 ≤ a i := Nat.succ_le_of_lt hgt
          -- from y ≥ 1, get y + 1 ≤ y + y
          have hy_le_2y : y + 1 ≤ y + y := by
            exact Nat.add_le_add_left hy_pos y
          -- from a i ≥ 2, get y + y ≤ a i * y
          have h2 : y * 2 ≤ a i * y := by
            -- first y * 2 ≤ y * a i
            have h2' : y * 2 ≤ y * a i := Nat.mul_le_mul_left y ha2
            -- rewrite right-hand side
            simpa [Nat.mul_comm] using h2'
          have h2y_ay : y + y ≤ a i * y := by
            simpa [nat_mul_two] using h2
          have hy_le_ay : y + 1 ≤ a i * y := le_trans hy_le_2y h2y_ay
          -- adding b i on the right preserves ≤
          have h_add : a i * y ≤ a i * y + b i := Nat.le_add_right _ _
          exact le_trans hy_le_ay h_add
      -- combine inequalities
      have : k + 2 ≤ a i * y + b i := le_trans hy_ge' step_ge
      exact this

theorem seqA_length (r : ℕ) (a b : Fin r → ℕ) (k : ℕ) : (seqA r a b k).length = r ^ k := by
  induction k with
  | zero =>
      simp [seqA]
  | succ k ih =>
      -- We prove a small helper lemma locally for lists of indices.
      have hbind_general (L : List (Fin r)) :
          (L >>= fun i => (seqA r a b k).map (fun x => a i * x + b i)).length
            = L.length * (seqA r a b k).length := by
        induction L with
        | nil =>
            simp
        | cons i L ihL =>
            have hi : ((seqA r a b k).map fun x => a i * x + b i).length = (seqA r a b k).length := by
              simp
            have hrest :
              (L >>= fun i => (seqA r a b k).map (fun x => a i * x + b i)).length
                = L.length * (seqA r a b k).length := ihL
            -- `simp` knows the definition of `>>=` for lists as `List.bind`.
            simp [hi, hrest, Nat.mul_succ, Nat.succ_mul, Nat.add_comm, Nat.add_left_comm, Nat.add_assoc] at *
      -- Now apply the general lemma to `finRange r`.
      have hlen_bind :
        (List.finRange r >>= fun i => (seqA r a b k).map (fun x => a i * x + b i)).length
          = (List.finRange r).length * (seqA r a b k).length := by
        exact hbind_general (List.finRange r)
      have hlen_finRange : (List.finRange r).length = r := by
        simpa using List.finRange_length r
      -- put everything together
      simp [seqA, ih, hlen_bind, hlen_finRange, Nat.pow_succ, Nat.mul_comm, Nat.mul_left_comm, Nat.mul_assoc]


theorem seqA_pos (r : ℕ) (a b : Fin r → ℕ) (ha : ∀ i, a i ≥ 1) (k : ℕ) : ∀ x ∈ seqA r a b k, x > 0 := by
  induction k with
  | zero =>
    intro x hx
    simp [seqA] at hx
    subst hx
    exact Nat.zero_lt_one
  | succ k ih =>
    intro x hx
    simp only [seqA, List.bind_eq_flatMap, List.mem_flatMap, List.mem_map] at hx
    rcases hx with ⟨i, _, y, hy, rfl⟩
    have h_pos_y := ih y hy
    have h_ge_a := ha i
    have h_mul : a i * y ≥ 1 := by
      rw [← Nat.mul_one 1]
      apply Nat.mul_le_mul h_ge_a
      exact Nat.succ_le_of_lt h_pos_y
    apply Nat.lt_of_lt_of_le (Nat.lt_of_succ_le h_mul)
    apply Nat.le_add_right

def seqA_q (r : ℕ) (a b : Fin r → ℕ) (k : ℕ) : List ℚ := (seqA r a b k).map (fun x => (x : ℚ))

theorem sorted_pos_get_ge_index (S : List ℕ) (hsort : S.Sorted (· < ·)) (hpos : ∀ x ∈ S, 0 < x) :
  ∀ i : Fin S.length, (i : ℕ) + 1 ≤ S.get i := by
  cases S with
    | nil =>
      intro i
      exact i.elim0
    | cons a as =>
      have hmono : StrictMono (a :: as).get := List.Sorted.get_strictMono hsort
      intro i
      induction i using Fin.induction with
      | zero =>
        have h0 : (a :: as).get ⟨0, Nat.zero_lt_succ as.length⟩ ∈ (a :: as) := List.get_mem _ _
        have hpos0 : 0 < (a :: as).get ⟨0, Nat.zero_lt_succ as.length⟩ := hpos _ h0
        exact Nat.succ_le_of_lt hpos0
      | succ i IH =>
        have hlt : Fin.castSucc i < i.succ := Fin.castSucc_lt_succ i
        have hmono_val : (a :: as).get (Fin.castSucc i) < (a :: as).get i.succ := hmono hlt
        calc
          (i.succ : ℕ) + 1 = (i : ℕ) + 2 := by simp
          _ ≤ (a :: as).get (Fin.castSucc i) + 1 := Nat.add_le_add_right IH 1
          _ ≤ (a :: as).get i.succ := Nat.succ_le_of_lt hmono_val

lemma sum_list_map_eq_sum_fin_custom {α β} [AddCommMonoid β] (L : List α) (f : α → β) :
  (L.map f).sum = ∑ i : Fin L.length, f (L.get i) := by
  nth_rw 1 [← List.ofFn_get L]
  rw [List.map_ofFn, List.sum_ofFn]
  rfl

lemma sum_fin_eq_sum_range_custom {β} [AddCommMonoid β] (n : ℕ) (f : ℕ → β) :
  ∑ i : Fin n, f i = ∑ i ∈ Finset.range n, f i := by
  rw [Fin.sum_univ_eq_sum_range]

theorem distinct_sum_bound_finset (S : Finset ℕ) (hpos : ∀ x ∈ S, x > 0) : ∑ x ∈ S, (1 : ℚ) / (x : ℚ) ≤ ∑ i ∈ Finset.range S.card, (1 : ℚ) / ((i : ℚ) + 1) := by
  classical
  let L := S.toList.insertionSort (· ≤ ·)
  have hL_perm : L.Perm S.toList := List.perm_insertionSort (· ≤ ·) S.toList
  have hL_sorted : L.Sorted (· ≤ ·) := List.sorted_insertionSort (· ≤ ·) S.toList
  have hL_len : L.length = S.card := by rw [List.Perm.length_eq hL_perm, S.length_toList]

  have hL_nodup : L.Nodup := List.Perm.nodup_iff hL_perm |>.mpr S.nodup_toList
  have hL_sorted_lt : L.Sorted (· < ·) := List.Sorted.lt_of_le hL_sorted hL_nodup
  have hL_pos : ∀ x ∈ L, 0 < x := by
    intro x hx
    have : x ∈ S := by rw [← Finset.mem_toList, ← hL_perm.mem_iff]; exact hx
    exact hpos x this

  have h_idx_le : ∀ i : Fin L.length, (i : ℕ) + 1 ≤ L.get i :=
    sorted_pos_get_ge_index L hL_sorted_lt hL_pos

  let g (k : ℕ) : ℚ := 1 / ((k : ℚ) + 1)

  calc
    ∑ x ∈ S, (1 : ℚ) / x = ∑ x ∈ S.toList.toFinset, (1 : ℚ) / x := by simp
    _ = (S.toList.map (fun x : ℕ => (1 : ℚ) / x)).sum := by
      rw [List.sum_toFinset]
      exact S.nodup_toList
    _ = (L.map (fun x : ℕ => (1 : ℚ) / x)).sum := by
      refine List.Perm.sum_eq (List.Perm.map _ hL_perm.symm)
    _ = ∑ i : Fin L.length, (1 : ℚ) / (L.get i) := by
      rw [sum_list_map_eq_sum_fin_custom]
    _ ≤ ∑ i : Fin L.length, g i := by
      apply Finset.sum_le_sum
      intro i _
      dsimp [g]
      apply one_div_le_one_div_of_le
      · norm_cast; exact Nat.succ_pos i
      · norm_cast; exact h_idx_le i
    _ = ∑ i ∈ Finset.range L.length, g i := by
      rw [sum_fin_eq_sum_range_custom]
    _ = ∑ i ∈ Finset.range S.card, (1 : ℚ) / ((i : ℚ) + 1) := by
      rw [hL_len]

theorem exists_ge_of_nodup_pos_length_ge (L : List ℕ) (hN : L.Nodup) (hpos : ∀ x ∈ L, x > 0) (M : ℕ) (hM1 : 1 ≤ M) (hMlen : M ≤ L.length) : ∃ x ∈ L, x ≥ M := by
  -- Use mergeSort with decidable (<=) relation
  let S := L.mergeSort (· ≤ ·)
  -- The correct name for the permutation lemma seems to be `List.mergeSort_perm`
  have hperm : S.Perm L := List.mergeSort_perm L (· ≤ ·)
  have hlen : S.length = L.length := hperm.length_eq
  have hnodup : S.Nodup := hN.perm hperm.symm
  -- Use `List.sorted_mergeSort'` which takes fewer explicit arguments according to search results
  have hsorted_le : S.Sorted (· ≤ ·) := List.sorted_mergeSort' (· ≤ ·) L

  -- S is strictly sorted because it is sorted (<=) and Nodup
  have hsorted_lt : S.Sorted (· < ·) := by
    rw [List.Sorted, List.pairwise_iff_get] at hsorted_le ⊢
    intro i j hij
    specialize hsorted_le i j hij
    have hne : S.get i ≠ S.get j := by
      -- Instead of rewriting, apply the iff directly or use logic
      intro heq
      rw [List.Nodup.get_inj_iff hnodup] at heq
      -- heq becomes i = j, which contradicts hij
      exact (ne_of_lt hij) heq
    exact Nat.lt_of_le_of_ne hsorted_le hne

  have hpos_S : ∀ x ∈ S, 0 < x := by
    intro x hx
    rw [hperm.mem_iff] at hx
    exact hpos x hx

  have hM_pos : 0 < M := Nat.lt_of_lt_of_le Nat.zero_lt_one hM1
  have hidx_lt : M - 1 < S.length := by
    rw [hlen]
    apply Nat.lt_of_lt_of_le _ hMlen
    apply Nat.sub_lt hM_pos Nat.zero_lt_one

  let i : Fin S.length := ⟨M - 1, hidx_lt⟩

  -- Apply the axiom
  have hge := sorted_pos_get_ge_index S hsorted_lt hpos_S i

  have hM_eq : (i : ℕ) + 1 = M := by
    simp only [i]
    exact Nat.sub_add_cancel hM1

  rw [hM_eq] at hge

  refine ⟨S.get i, ?_, hge⟩
  rw [← hperm.mem_iff]
  exact List.get_mem _ _

theorem nodup_nat_lt_bound_length_le (M : ℕ) (L : List ℕ) (hN : L.Nodup) (hLt : ∀ x ∈ L, x < M) : L.length ≤ M := by
  let L1 := L.map (· + 1)
  have hN1 : L1.Nodup := List.Nodup.map (fun x y h => Nat.succ.inj h) hN
  let S := L1.mergeSort (· ≤ ·)
  have h_perm : List.Perm S L1 := List.mergeSort_perm L1 (· ≤ ·)
  have hS_len : S.length = L.length := by rw [h_perm.length_eq, List.length_map]

  -- Use the Mathlib wrapper for sortedness
  have hS_sorted : S.Sorted (· ≤ ·) := List.sorted_mergeSort' (· ≤ ·) L1

  have hS_nodup : S.Nodup := h_perm.nodup_iff.mpr hN1
  have hS_strict : S.Sorted (· < ·) := List.Sorted.lt_of_le hS_sorted hS_nodup
  have hS_pos : ∀ x ∈ S, 0 < x := by
    intro x hx
    rw [h_perm.mem_iff, List.mem_map] at hx
    rcases hx with ⟨y, _, rfl⟩
    exact Nat.succ_pos y

  cases h_S_nil : S with
  | nil =>
    rw [h_S_nil, List.length_nil] at hS_len
    rw [← hS_len]
    exact Nat.zero_le M
  | cons head tail =>
    let last_idx : Fin S.length := ⟨S.length - 1, Nat.sub_lt (by
      rw [h_S_nil, List.length_cons]
      apply Nat.zero_lt_succ) (by simp)⟩

    have h_axiom := sorted_pos_get_ge_index S hS_strict hS_pos last_idx

    let last := S.get last_idx
    have h_last_mem : last ∈ S := List.get_mem S last_idx

    have h_last_lt : last < M + 1 := by
      rw [h_perm.mem_iff, List.mem_map] at h_last_mem
      rcases h_last_mem with ⟨y, hyL, heq⟩
      rw [← heq]
      apply Nat.succ_lt_succ
      exact hLt y hyL

    have h_idx_val : (last_idx : ℕ) = L.length - 1 := by
      simp only [last_idx]
      rw [hS_len]

    rw [h_idx_val] at h_axiom

    have h_len_pos : L.length > 0 := by
      rw [← hS_len, h_S_nil]
      exact Nat.zero_lt_succ _

    have : L.length - 1 + 1 = L.length := Nat.sub_add_cancel h_len_pos
    rw [this] at h_axiom

    -- h_axiom : L.length ≤ last
    apply Nat.le_trans h_axiom
    exact Nat.le_of_lt_succ h_last_lt

theorem sum_filter_map_eq_sum_map_ite {α β : Type} [AddCommMonoid β] (L : List α) (p : α → Prop) [DecidablePred p] (f : α → β) :
  ((L.filter p).map f).sum = (L.map (fun x => if p x then f x else 0)).sum := by
  induction L with
  | nil => simp
  | cons head tail ih =>
    simp only [List.map_cons, List.sum_cons]
    by_cases h : p head
    · rw [if_pos h]
      rw [List.filter_cons]
      simp [h, ih]
    · rw [if_neg h]
      rw [List.filter_cons]
      simp [h, ih]

def sum_inv (L : List ℕ) : ℚ := (L.map (fun x => (1 : ℚ) / (x : ℚ))).sum

lemma list_sum_map_eq_finset_sum {α β : Type*} [DecidableEq α] [AddCommMonoid β]
  (l : List α) (f : α → β) (h : l.Nodup) :
  (l.map f).sum = ∑ x ∈ l.toFinset, f x := by
  rw [← List.sum_toFinset _ h]

theorem distinct_sum_bound (L : List ℕ) (h : L.Nodup) (hpos : ∀ x ∈ L, x > 0) :
  (L.map (fun x => (1 : ℚ) / (x : ℚ))).sum ≤ ((List.range L.length).map (fun x => (1 : ℚ) / ((x + 1) : ℚ))).sum := by
  -- Main proof
  -- Fix `List.nodup_range` (again). It is NOT a function, it's a theorem `(List.range n).Nodup`.
  -- `List.nodup_range L.length` is valid ONLY if `List.nodup_range` takes an argument.
  -- In Mathlib 4, it is `List.nodup_range (n : Nat) : (List.range n).Nodup`.
  -- Wait, the error says: `Function expected at List.nodup_range but this term has type (List.range ?m).Nodup`.
  -- This means `List.nodup_range` is a term (proof), not a function.
  -- It means the type is ALREADY specialized or inferred?
  -- Or `List.nodup_range` refers to the proof for implicit arguments?
  -- Ah, `List.nodup_range` is `forall n, ...` or just `(List.range n).Nodup`?
  -- If it's a global theorem `forall n`, I apply it.
  -- If it's a field or something...

  -- Actually, `List.nodup_range` usually takes `n`.
  -- Let's try omitting the argument and let Lean infer it from context?
  -- Or explicit `@List.nodup_range L.length`.

  -- The error `List.nodup_range` having type `(List.range ?m).Nodup` suggests it takes NO explicit arguments and relies on implicit `n`.
  -- So I should just write `List.nodup_range`.

  -- Regarding `convert` failure:
  -- The goals generated are:
  -- 1. `L.length = S.card` (from matching indices in sum)
  -- 2. `∀ x ∈ S, x > 0` (argument to distinct_sum_bound_finset)

  -- I need to be careful with `convert`.
  -- `convert distinct_sum_bound_finset S _ using 1` might be safer.

  have h_flatmap_eq_map {α β : Type} (l : List α) (f : α → β) :
    (l.flatMap (fun a => [f a])) = l.map f := by
    induction l with
    | nil => rfl
    | cons head tail ih => simp [ih]

  -- Step 1: LHS
  let S := L.toFinset
  let sum_S := ∑ x ∈ S, (1 : ℚ) / (x : ℚ)
  let G_lhs := (L.map (fun x => (1 : ℚ) / (x : ℚ))).sum
  let f_comp := (fun x => (1 : ℚ) / x) ∘ (fun (a : ℕ) => (a : ℚ))

  have h_lhs : G_lhs = sum_S := by
    dsimp [G_lhs]
    rw [h_flatmap_eq_map L (fun a => (a : ℚ))]
    rw [List.map_map]
    rw [list_sum_map_eq_finset_sum L f_comp h]
    rfl

  -- Step 2: RHS
  let R := List.range L.length
  let sum_R := ∑ x ∈ Finset.range L.length, (1 : ℚ) / ((x : ℚ) + 1)
  let G_rhs := (R.map (fun x => (1 : ℚ) / ((x + 1) : ℚ))).sum
  let g_comp := (fun x => (1 : ℚ) / (x + 1)) ∘ (fun (a : ℕ) => (a : ℚ))

  -- Try without argument for nodup_range
  have hR : R.Nodup := List.nodup_range

  have h_rhs : G_rhs = sum_R := by
    dsimp [G_rhs]
    rw [h_flatmap_eq_map R (fun a => (a : ℚ))]
    rw [List.map_map]
    rw [list_sum_map_eq_finset_sum R g_comp hR]
    rw [List.toFinset_range]
    rfl

  change G_lhs ≤ G_rhs
  rw [h_lhs, h_rhs]

  have h_card : S.card = L.length := List.toFinset_card_of_nodup h

  -- Explicitly construct the proof term to avoid `convert` ambiguity
  have h_axiom := distinct_sum_bound_finset S (by
    intro x hx
    apply hpos
    exact List.mem_toFinset.mp hx
  )

  -- h_axiom : sum_S ≤ ∑ i ∈ range S.card, ...
  rw [h_card] at h_axiom
  exact h_axiom

noncomputable def sum_inv_ge (M : ℕ) (L : List ℕ) : ℝ :=
  ((L.filter (fun x => M ≤ x)).map (fun x => (1 : ℝ) / (x : ℝ))).sum

lemma flatMap_pure_eq_map' {α β : Type} (l : List α) (f : α → β) : l.flatMap (fun x => [f x]) = l.map f := by induction l <;> simp [*]

theorem sum_inv_ge_eq_sum_ite (M : ℕ) (L : List ℕ) :
  sum_inv_ge M L =
    (L.map (fun x => if M ≤ x then (1 : ℝ) / (x : ℝ) else 0)).sum := by
  dsimp [sum_inv_ge]
  simp only [flatMap_pure_eq_map']
  rw [List.map_map]
  exact sum_filter_map_eq_sum_map_ite L (fun x => M ≤ x) (fun x => (1 : ℝ) / (x : ℝ))

lemma list_sum_pos_of_forall_pos :
  ∀ (L : List ℝ), (∀ x ∈ L, 0 < x) → L ≠ [] → 0 < L.sum
  | [], hpos, hne => by
      cases hne rfl
  | a :: L, hpos, hne => by
      have hapos : 0 < a := hpos a (by simp)
      by_cases hLnil : L = []
      · subst hLnil
        simpa [List.sum_cons] using hapos
      · have hpos_tail : ∀ x ∈ L, 0 < x := by
          intro x hx
          apply hpos x
          exact List.mem_cons_of_mem _ hx
        have ih : 0 < L.sum := list_sum_pos_of_forall_pos L hpos_tail hLnil
        have : 0 < a + L.sum := add_pos hapos ih
        simpa [List.sum_cons] using this

theorem sum_inv_ge_pos_of_exists_ge (M : ℕ) (L : List ℕ) (hpos : ∀ x ∈ L, x > 0) (hex : ∃ x ∈ L, M ≤ x) : sum_inv_ge M L > 0 := by
  unfold sum_inv_ge
  let S := L.filter (fun x => M ≤ x)

  apply list_sum_pos_of_forall_pos
  · intro y hy
    rw [List.mem_map] at hy
    rcases hy with ⟨x_real, hx_real, rfl⟩
    simp at hx_real
    rcases hx_real with ⟨n, hnS, hnx⟩
    rw [hnx]

    have hn_pos : 0 < n := by
      exact hpos n hnS.1

    rw [one_div_pos]
    exact_mod_cast hn_pos

  · intro h_nil
    simp only [List.map_eq_nil_iff] at h_nil

    rcases hex with ⟨x0, hx0L, hx0M⟩
    have hx0S : x0 ∈ S := by
      rw [List.mem_filter]
      simp [hx0L, hx0M]

    have h_mem : (x0 : ℝ) ∈ (do let a ← S; pure (a : ℝ)) := by
      simp
      -- The error `Invalid ⟨...⟩ notation: expected List.Mem x0 S has more than one constructor`
      -- implies the goal WAS `x0 ∈ S` (List.Mem), NOT `Exists ...`.
      -- This means `simp` simplified the existential away?
      -- `∃ a, a ∈ S ∧ ↑a = ↑x0`.
      -- If `simp` knows `↑a = ↑x0 <-> a = x0` (since coe is injective for Nat->Real),
      -- then `simp` simplifies to `∃ a, a ∈ S ∧ a = x0`, which is `x0 ∈ S`.
      -- So the goal became `x0 ∈ S`.
      -- And `refine ⟨x0, ...⟩` was trying to build a proof of `x0 ∈ S` (which is `List.Mem`),
      -- but `List.Mem` has `head` and `tail` constructors, so `<>` syntax fails unless it's a structure.
      -- So we just need to prove `x0 ∈ S`.
      exact hx0S

    rw [h_nil] at h_mem
    contradiction

theorem seqA_large_part_eventually_positive (r : ℕ) (a b : Fin r → ℕ) (ha : ∀ i, a i ≥ 1) (hr : 2 ≤ r) (h_nodup_all : ∀ k, (seqA r a b k).Nodup) (M0 : ℕ) :
  ∃ K : ℕ, sum_inv_ge M0 (seqA r a b K) > 0 := by
  classical
  cases M0 with
  | zero =>
    -- Case M0 = 0
    refine ⟨0, ?_⟩
    -- All elements in seqA r a b 0 are positive
    have hpos0 : ∀ x ∈ seqA r a b 0, x > 0 := by
      intro x hx
      exact seqA_pos r a b ha 0 x hx
    -- There exists an element ≥ 0 (take x = 1 in [1])
    have hex0 : ∃ x ∈ seqA r a b 0, 0 ≤ x := by
      refine ⟨1, ?_, ?_⟩
      · simp [seqA]
      · norm_num
    -- Apply positivity lemma
    have := sum_inv_ge_pos_of_exists_ge 0 (seqA r a b 0) hpos0 hex0
    simpa using this
  | succ m =>
    -- Case M0 = m+1 ≥ 1
    -- Introduce the helper lemma locally (already available in context)
    have h_two_pow_ge : m.succ ≤ 2 ^ m.succ := by
      -- Prove m+1 ≤ 2*(m+1) ≤ 2^(m+1)
      have h1 : m.succ ≤ 2 * m.succ := by
        have h12 : (1 : ℕ) ≤ 2 := by decide
        simpa [one_mul, two_mul] using Nat.mul_le_mul_right m.succ h12
      have h2 : 2 * m.succ ≤ 2 ^ m.succ := by
        -- Note: m.succ = (m+1), so this matches Nat.mul_le_pow
        have := Nat.mul_le_pow (by decide : (2 : ℕ) ≠ 1) m.succ
        simpa [Nat.succ_eq_add_one, two_mul] using this
      exact le_trans h1 h2
    -- Choose K = m+1
    let K : ℕ := m.succ
    let L : List ℕ := seqA r a b K
    have hN : L.Nodup := by
      simpa [L, K] using h_nodup_all K
    have hpos : ∀ x ∈ L, x > 0 := by
      intro x hx
      have hx' : x ∈ seqA r a b K := by simpa [L] using hx
      exact seqA_pos r a b ha K x hx'
    -- Show length of L is at least m+1 using growth of powers
    have h_pow_ge : m.succ ≤ r ^ m.succ := by
      -- First, m+1 ≤ 2^(m+1)
      have h1 : m.succ ≤ 2 ^ m.succ := h_two_pow_ge
      -- Then, 2^(m+1) ≤ r^(m+1) since r ≥ 2
      have h2 : 2 ^ m.succ ≤ r ^ m.succ := by
        have hbase : (2 : ℕ) ≤ r := hr
        -- monotonicity of powers in the base
        have hpow := pow_le_pow_left' (M := ℕ) hbase
        exact hpow m.succ
      exact le_trans h1 h2
    have hlen : m.succ ≤ L.length := by
      -- Use seqA_length to relate length to r^(m+1)
      have : m.succ ≤ r ^ m.succ := h_pow_ge
      simpa [L, K, seqA_length] using this
    have hM1 : 1 ≤ m.succ := Nat.succ_le_succ (Nat.zero_le m)
    -- Use combinatorial lemma to find an element ≥ m+1 in L
    have hex_int : ∃ x ∈ L, x ≥ m.succ := by
      have := exists_ge_of_nodup_pos_length_ge L hN hpos m.succ hM1 hlen
      simpa using this
    -- Now apply positivity lemma for sum_inv_ge
    have hsum_pos := sum_inv_ge_pos_of_exists_ge m.succ L hpos hex_int
    refine ⟨K, ?_⟩
    simpa [K, L] using hsum_pos

def sum_inv_q (L : List ℚ) : ℚ := (L.map (fun x => 1 / x)).sum

theorem sum_inv_real_eq (L : List ℕ) : (((sum_inv L : ℚ)) : ℝ) = (L.map (fun x => (1 : ℝ) / (x : ℝ))).sum := by
  -- First prove a helper lemma describing `sum_inv` on a cons list.
  have sum_inv_cons : ∀ (x : ℕ) (xs : List ℕ),
      sum_inv (x :: xs) = (1 : ℚ) / (x : ℚ) + sum_inv xs := by
    intro x xs
    simp [sum_inv, List.map_cons, List.sum_cons]
  -- Now prove the main statement by induction on `L`.
  induction L with
  | nil =>
    simp [sum_inv]
  | cons x xs ih =>
    -- Expand `sum_inv` on the left using the helper lemma.
    have hL : sum_inv (x :: xs) = (1 : ℚ) / (x : ℚ) + sum_inv xs := sum_inv_cons x xs
    -- Rewrite the left-hand side.
    -- Original goal:
    --   ↑(sum_inv (x :: xs)) = ((x :: xs).map (fun x => (1 : ℝ) / (x : ℝ))).sum
    -- After rewriting:
    --   ↑((1 : ℚ) / (x : ℚ) + sum_inv xs)
    --     = (1 : ℝ) / (x : ℝ) + (xs.map (fun x => (1 : ℝ) / (x : ℝ))).sum
    -- We perform this rewrite and simplify.
    have := hL
    -- Now change the goal accordingly.
    -- Use `simp` with `Rat.cast_add` and the induction hypothesis.
    -- Apply `simp` directly to the goal.
    simp [hL, sum_inv_cons, ih, Rat.cast_add, one_div]


theorem sum_inv_ge_le_sum_inv (M : ℕ) (L : List ℕ) (hpos : ∀ x ∈ L, x > 0) :
  sum_inv_ge M L ≤ (((sum_inv L : ℚ)) : ℝ) := by
  classical
  -- Reduce the right-hand side to the real sum over the whole list
  suffices hmain :
      sum_inv_ge M L ≤ (L.map (fun x => (1 : ℝ) / (x : ℝ))).sum by
    simpa [sum_inv_real_eq] using hmain
  -- Prove the inequality by induction on L, carrying the positivity hypothesis
  revert hpos
  induction L with
  | nil =>
      intro hpos
      simp [sum_inv_ge]
  | cons x xs ih =>
      intro hpos
      -- positivity of head and tail elements
      have hx_pos_nat : x > 0 := hpos x (by simp)
      have hxs_pos : ∀ y ∈ xs, y > 0 := by
        intro y hy
        exact hpos y (by simp [hy])
      have ih' : sum_inv_ge M xs ≤ (xs.map (fun t => (1 : ℝ) / (t : ℝ))).sum :=
        ih hxs_pos
      -- positivity of 1 / x in ℝ
      have hx_pos_real : (0 : ℝ) < (x : ℝ) := by
        exact_mod_cast hx_pos_nat
      have hx_inv_nonneg : 0 ≤ (1 : ℝ) / (x : ℝ) :=
        le_of_lt (one_div_pos.mpr hx_pos_real)
      -- case split on whether x survives the filter
      by_cases hMx : M ≤ x
      · -- head is kept by the filter, add the same positive term to both sides
        have hadd :
            (1 : ℝ) / (x : ℝ) + sum_inv_ge M xs ≤
              (1 : ℝ) / (x : ℝ) + (xs.map (fun t => (1 : ℝ) / (t : ℝ))).sum :=
          add_le_add_left ih' _
        simpa [sum_inv_ge, hMx] using hadd
      · -- head is dropped by the filter, the full sum gains an extra positive term
        have hsum_le :
            (xs.map (fun t => (1 : ℝ) / (t : ℝ))).sum ≤
              (1 : ℝ) / (x : ℝ) + (xs.map (fun t => (1 : ℝ) / (t : ℝ))).sum :=
          le_add_of_nonneg_left hx_inv_nonneg
        have htrans :
            sum_inv_ge M xs ≤
              (1 : ℝ) / (x : ℝ) + (xs.map (fun t => (1 : ℝ) / (t : ℝ))).sum :=
          le_trans ih' hsum_le
        simpa [sum_inv_ge, hMx] using htrans

theorem sum_inv_recurrence (r : ℕ) (a b : Fin r → ℕ) (A : List ℕ) :
  ((List.finRange r >>= (fun i => A.map (fun x => a i * x + b i))).map (fun x => (1 : ℚ) / (x : ℚ))).sum = (A.map (fun x => (∑ i : Fin r, (1 : ℚ) / ((a i * x + b i) : ℚ)))).sum := by
  classical
  simp [List.bind_eq_flatMap, list_bind_sum_map, flatMap_sum_map,
        list_sum_finRange_eq_finset_sum, list_sum_finset_sum_comm,
        List.map_map, Function.comp]
  refine Finset.sum_congr rfl ?_
  intro i hi
  apply congrArg (fun L : List ℚ => L.sum)
  -- prove equality of maps by rewriting with `List.map_eq_map` via `List.map_eq_pmap`?
  -- Instead, we can just use `rfl` because `simp` already normalized both sides.
  -- Check if the goal is already solved by `simp`.
  simp [Function.comp, Nat.cast_add, Nat.cast_mul, mul_comm, mul_left_comm,
        mul_assoc, add_comm, add_left_comm, add_assoc]

theorem sum_inv_recurrence_pointwise (r : ℕ) (a b : Fin r → ℕ) (ha : ∀ i, a i ≥ 1) (x : ℕ) (hxpos : x > 0) :
  (∑ i : Fin r, (1 : ℚ) / ((a i * x + b i : ℕ) : ℚ))
    ≥ (∑ i : Fin r, (1 : ℚ) / (a i : ℚ)) * ((1 : ℚ) / (x : ℚ))
      - (∑ i : Fin r, (b i : ℚ) / (a i : ℚ)^2) * ((1 : ℚ) / (x : ℚ)^2) := by
  have hx_q : (x : ℚ) > 0 := Nat.cast_pos.mpr hxpos
  rw [Finset.sum_mul, Finset.sum_mul, ← Finset.sum_sub_distrib]
  apply Finset.sum_le_sum
  intro i _
  have ha_q : (a i : ℚ) > 0 := Nat.cast_pos.mpr (Nat.lt_of_succ_le (ha i))
  have hb_q : (b i : ℚ) ≥ 0 := Nat.cast_nonneg (b i)
  have h := inv_add_ge_sub_square_helper (a i) x (b i) hx_q ha_q hb_q
  simp only [ge_iff_le] at h ⊢
  convert h using 1
  · field_simp
  · push_cast
    rfl

theorem per_point_recurrence_large_x (r : ℕ) (a b : Fin r → ℕ) (ha : ∀ i, a i ≥ 1) (hsum : ∑ i : Fin r, (1 : ℚ) / (a i : ℚ) > 1) : ∃ (ρ : ℝ) (hρ : ρ > 1) (M0 : ℕ), 1 ≤ M0 ∧ ∀ n : ℕ, n ≥ M0 → (∑ i : Fin r, (1 : ℝ) / ((a i * n + b i : ℕ) : ℝ)) ≥ ρ * ((1 : ℝ) / (n : ℝ)) := by
  -- Define rational sums
  let σ_Q := ∑ i : Fin r, (1 : ℚ) / (a i : ℚ)
  let C_Q := ∑ i : Fin r, (b i : ℚ) / (a i : ℚ)^2

  -- Real versions
  let σ_R : ℝ := σ_Q
  let C_R : ℝ := C_Q

  -- σ_R > 1
  have hσ_R : σ_R > (1 : ℝ) := by
    -- cast the rational inequality hsum to ℝ
    have hσ_Q : (σ_Q : ℝ) > (1 : ℝ) := by
      exact_mod_cast hsum
    simpa [σ_R] using hσ_Q

  -- Choose ρ with 1 < ρ < σ_R
  obtain ⟨ρ, hρ_gt_1, hρ_lt_σ⟩ := exists_between hσ_R

  -- Define δ and show positivity
  let δ : ℝ := σ_R - ρ
  have hδ : δ > 0 := sub_pos.mpr hρ_lt_σ

  -- Use Archimedean property to bound C_R by n * δ for large n
  obtain ⟨N, hN⟩ : ∃ N : ℕ, C_R / δ ≤ N := exists_nat_ge (C_R / δ)

  -- Define threshold M0
  let M0 : ℕ := max N 1

  refine ⟨ρ, hρ_gt_1, M0, ?_, ?_⟩
  · -- M0 ≥ 1
    have : (1 : ℕ) ≤ M0 := by
      unfold M0
      exact le_max_right _ _
    exact this

  · -- Main inequality for all n ≥ M0
    intro n hn

    -- n ≥ 1, hence n > 0 as Nat and as ℝ
    have hM0_ge1 : (1 : ℕ) ≤ M0 := by
      unfold M0
      exact le_max_right _ _
    have hn_ge1 : (1 : ℕ) ≤ n := le_trans hM0_ge1 hn
    have hn_pos_nat : 0 < n := Nat.succ_le_iff.mp hn_ge1
    have hn_pos_real : (0 : ℝ) < (n : ℝ) := by exact_mod_cast hn_pos_nat

    -- Bound condition: (n : ℝ) * (σ_R - ρ) ≥ C_R
    have h_bound_condition : (n : ℝ) * (σ_R - ρ) ≥ C_R := by
      -- First, N ≤ n as naturals, hence as reals
      have hN_le_n_nat : N ≤ n := by
        -- from max N 1 ≤ n get N ≤ n
        have hM0_le_n : M0 ≤ n := hn
        have hN_le_M0 : N ≤ M0 := by
          unfold M0
          exact le_max_left _ _
        exact le_trans hN_le_M0 hM0_le_n
      have hN_le_n_real : (N : ℝ) ≤ (n : ℝ) := by exact_mod_cast hN_le_n_nat

      -- From hN : C_R / δ ≤ N (viewed in ℝ)
      have hN_real : C_R / δ ≤ (N : ℝ) := hN

      -- Multiply inequalities by δ > 0
      have h1 : C_R / δ * δ ≤ (N : ℝ) * δ :=
        mul_le_mul_of_nonneg_right hN_real (le_of_lt hδ)
      have h2 : (N : ℝ) * δ ≤ (n : ℝ) * δ :=
        mul_le_mul_of_nonneg_right hN_le_n_real (le_of_lt hδ)

      -- Rewrite C_R / δ * δ = C_R
      have hδ_ne : (δ : ℝ) ≠ 0 := ne_of_gt hδ
      have hC_eq : C_R / δ * δ = C_R := by
        field_simp [δ, hδ_ne]

      -- Chain inequalities to get C_R ≤ (n : ℝ) * δ
      have : C_R ≤ (n : ℝ) * δ := by
        calc
          C_R = C_R / δ * δ := by simpa [hC_eq]
          _ ≤ (N : ℝ) * δ := h1
          _ ≤ (n : ℝ) * δ := h2
      -- Replace δ by σ_R - ρ
      simpa [δ, mul_comm, mul_left_comm, mul_assoc] using this

    -- Now convert bound condition to an inequality of fractions
    have h_real_imp : σ_R / (n : ℝ) - C_R / (n : ℝ)^2 ≥ ρ / (n : ℝ) := by
      -- Use the axiom div_sub_div_sq_ge_iff_mul
      have := (div_sub_div_sq_ge_iff_mul σ_R ρ C_R (n : ℝ) hn_pos_real).mpr h_bound_condition
      -- Just change notation n:ℝ vs (n:ℝ)
      simpa using this

    -- Rational inequality from axiom, at integer n
    have h_rat_ineq := sum_inv_recurrence_pointwise r a b ha n hn_pos_nat

    -- Work in reals: define the real sum S_n
    let S_n : ℝ := ∑ i : Fin r, (1 : ℝ) / ((a i * n + b i : ℕ) : ℝ)

    -- Show S_n ≥ σ_R / n - C_R / n^2
    have h_S_ge : S_n ≥ σ_R / (n : ℝ) - C_R / (n : ℝ)^2 := by
      -- It is more convenient to work with the ≤ form
      change σ_R / (n : ℝ) - C_R / (n : ℝ)^2 ≤ S_n

      -- Cast the rational inequality to ℝ
      have h_rat_ineq_R :
          ((∑ i : Fin r, (1 : ℚ) / ((a i * n + b i : ℕ) : ℚ)) : ℝ) ≥
            ((∑ i : Fin r, (1 : ℚ) / (a i : ℚ)) * (1 / (n : ℚ)) -
              (∑ i : Fin r, (b i : ℚ) / (a i : ℚ)^2) * (1 / (n : ℚ)^2) : ℝ) := by
        exact_mod_cast h_rat_ineq

      -- Reinterpret it as a ≤ inequality
      have h_rat_ineq_le :
          ((∑ i : Fin r, (1 : ℚ) / (a i : ℚ)) * (1 / (n : ℚ)) -
              (∑ i : Fin r, (b i : ℚ) / (a i : ℚ)^2) * (1 / (n : ℚ)^2) : ℝ) ≤
            ((∑ i : Fin r, (1 : ℚ) / ((a i * n + b i : ℕ) : ℚ)) : ℝ) := by
        simpa using h_rat_ineq_R

      -- From a - b ≤ c, deduce a ≤ c + b
      have h_alt :
          ((∑ i : Fin r, (1 : ℚ) / (a i : ℚ)) * (1 / (n : ℚ)) : ℝ) ≤
            ((∑ i : Fin r, (1 : ℚ) / ((a i * n + b i : ℕ) : ℚ)) : ℝ) +
              ((∑ i : Fin r, (b i : ℚ) / (a i : ℚ)^2) * (1 / (n : ℚ)^2) : ℝ) := by
        -- apply sub_le_iff_le_add to h_rat_ineq_le
        have := (sub_le_iff_le_add).mp h_rat_ineq_le
        -- Now `this` has the desired type
        simpa [add_comm, add_left_comm, add_assoc] using this

      -- Now rewrite everything in terms of S_n, σ_R, C_R and ℝ-division by n
      have h_alt' : σ_R / (n : ℝ) ≤ S_n + C_R / (n : ℝ)^2 := by
        -- simplify both sides
        simpa [S_n, σ_R, C_R, σ_Q, C_Q, one_div, div_eq_mul_inv, pow_two] using h_alt

      -- Finally, go back from a ≤ b + c to a - c ≤ b
      have h_sub : σ_R / (n : ℝ) - C_R / (n : ℝ)^2 ≤ S_n :=
        (sub_le_iff_le_add).2 h_alt'
      exact h_sub

    -- Finally combine inequalities and rewrite the goal
    have h_final : S_n ≥ ρ * (1 / (n : ℝ)) := by
      -- From S_n ≥ σ_R/n - C_R/n^2 and σ_R/n - C_R/n^2 ≥ ρ/n
      have h_chain : S_n ≥ ρ / (n : ℝ) := by
        -- S_n ≥ A and A ≥ B implies S_n ≥ B
        exact ge_trans h_S_ge h_real_imp
      -- Rewrite ρ / n as ρ * (1 / n)
      simpa [one_div, div_eq_mul_inv] using h_chain

    -- The goal is about the same S_n
    simpa [S_n, one_div, div_eq_mul_inv] using h_final

lemma algebraic_helper (y b : ℚ) (hy : y > 0) (hb : b ≥ 0) :
  1 / (y + b) ≥ 1 / y - b / (y ^ 2) := by
  have y_sq_pos : 0 < y ^ 2 := pow_pos hy 2
  have y_plus_b_pos : 0 < y + b := add_pos_of_pos_of_nonneg hy hb
  rw [ge_iff_le, sub_le_iff_le_add]
  field_simp
  nlinarith

lemma map_flatMap {α β γ} (l : List α) (f : α → List β) (g : β → γ) :
  (l.flatMap f).map g = l.flatMap (fun x => (f x).map g) := by
  induction l <;> simp [*]

lemma sum_flatMap {r : ℕ} {α} [AddMonoid α] (l : List (Fin r)) (f : Fin r → List α) :
  (l.flatMap f).sum = (l.map (fun x => (f x).sum)).sum := by
  induction l <;> simp [*, List.sum_append]

lemma sum_map_finRange {r : ℕ} {α} [AddCommMonoid α] (f : Fin r → α) :
  ((List.finRange r).map f).sum = ∑ i : Fin r, f i := by
  simp [Finset.sum, Finset.univ, Fintype.elems, Fin.fintype]

lemma sum_sub_map {α} [AddCommGroup α] (L : List ℕ) (f g : ℕ → α) :
  (L.map (fun x => f x - g x)).sum = (L.map f).sum - (L.map g).sum := by
  induction L <;> simp [*, add_sub_add_comm]

lemma list_sum_mul (L : List ℕ) (g : ℕ → ℚ) (c : ℚ) :
  c * (L.map g).sum = (L.map (fun x => c * g x)).sum := by
  induction L <;> simp [*, mul_add]

lemma sum_le_sum_map (L : List ℕ) (f g : ℕ → ℚ) (h : ∀ x ∈ L, f x ≥ g x) :
  (L.map f).sum ≥ (L.map g).sum := by
  induction L
  case nil => simp
  case cons head tail ih =>
    simp
    apply add_le_add
    · apply h; simp
    · apply ih; intro x hx; apply h; simp [hx]

lemma map_map_special {α β γ} (L : List α) (f : α → β) (g : β → γ) :
  (L.map f).map g = L.map (fun x => g (f x)) := by
  rw [List.map_map]; rfl

lemma flatMap_pure_eq_map'' {α β : Type} (l : List α) (f : α → β) :
  l.flatMap (fun x => [f x]) = l.map f := by
  induction l <;> simp [*]

theorem sum_inv_recurrence_bound (r : ℕ) (a b : Fin r → ℕ) (ha : ∀ i, a i ≥ 1) (A : List ℕ) (h_nodup : A.Nodup) (hpos : ∀ x ∈ A, x > 0) : ((List.finRange r >>= (fun i => A.map (fun x => a i * x + b i))).map (fun x => (1 : ℚ) / (x : ℚ))).sum ≥ (∑ i : Fin r, (1 : ℚ) / (a i : ℚ)) * ((A.map (fun x => (1 : ℚ) / (x : ℚ))).sum) - (∑ i : Fin r, (b i : ℚ) / (a i : ℚ)^2) * ((A.map (fun x => (1 : ℚ) / (x : ℚ)^2)).sum) := by
  -- Start main proof
  rw [sum_inv_recurrence r a b A]
  simp only [List.pure_def, List.bind_eq_flatMap, flatMap_pure_eq_map]
  simp only [List.map_map, Function.comp]
  rw [list_sum_mul, list_sum_mul]
  rw [←sum_sub_map]
  apply sum_le_sum_map
  intro x hx

  -- Pointwise inequality
  dsimp only [Function.comp, Nat.cast]
  rw [div_eq_mul_inv _ (x:ℚ), Finset.sum_mul]
  rw [div_eq_mul_inv _ ((x:ℚ)^2), Finset.sum_mul]
  rw [←Finset.sum_sub_distrib]

  apply Finset.sum_le_sum
  intro i hi

  have x_pos : (x : ℚ) > 0 := by norm_cast; apply hpos _ hx
  have ai_pos : (a i : ℚ) > 0 := by norm_cast; linarith [ha i]
  have y_pos : (a i : ℚ) * x > 0 := mul_pos ai_pos x_pos
  have bi_nonneg : (b i :ℚ) ≥ 0 := by norm_cast; apply zero_le

  have h_alg := algebraic_helper ((a i : ℚ) * x) (b i) y_pos bi_nonneg

  -- Normalize everything and close
  simp at h_alg ⊢
  rw [pow_two] at h_alg ⊢
  field_simp at h_alg ⊢
  convert h_alg using 1

def sum_inv_sq (L : List ℕ) : ℚ := (L.map (fun x => (1 : ℚ) / (x : ℚ)^2)).sum

theorem sum_inv_sq_le_two (L : List ℕ) (h_nodup : L.Nodup) (hpos : ∀ x ∈ L, x > 0) :
  (L.map (fun x => (1 : ℚ) / (x : ℚ)^2)).sum ≤ (2 : ℚ) := by
  have h_eq : (List.map (fun x => (1 : ℚ) / x^2) (do let a <- L; pure ↑a)) = List.map (fun n : ℕ => 1 / (n : ℚ)^2) L := by
    rw [List.bind_eq_flatMap]
    generalize L = l
    induction l with
    | nil => rfl
    | cons head tail ih =>
      simp only [List.pure_def] at ih ⊢
      simp only [List.flatMap_cons, List.map_cons, List.map_append]
      simp only [List.map_nil, List.singleton_append]
      rw [ih]

  rw [h_eq]

  let f (n : ℕ) : ℚ := 1 / (n : ℚ)^2
  let S := L.toFinset

  have h_bound : ∑ x ∈ S, f x ≤ 2 := by
    by_cases h : S = ∅
    · rw [h, Finset.sum_empty]
      norm_num
    · have hne : S.Nonempty := Finset.nonempty_of_ne_empty h
      let m := S.max' hne
      have h_sub : S ⊆ Finset.Ioo 0 (m + 1) := by
        intro x hx
        rw [Finset.mem_Ioo]
        constructor
        · rw [List.mem_toFinset] at hx
          exact hpos x hx
        · exact Nat.lt_succ_of_le (S.le_max' x hx)
      refine le_trans (Finset.sum_le_sum_of_subset_of_nonneg h_sub ?_) ?_
      · intro x _ _
        dsimp [f]
        positivity
      · have bound := sum_Ioo_inv_sq_le (α := ℚ) 0 (m + 1)
        simp only [Nat.cast_zero, zero_add, div_one] at bound
        convert bound using 1
        congr
        ext x
        dsimp [f]
        simp only [inv_eq_one_div]
        -- norm_cast was redundant here, removing it completely.

  rw [List.sum_toFinset f h_nodup] at h_bound
  exact h_bound

theorem sum_map_ite_bound {α : Type} (L : List α) (f g : α → ℝ) (p : α → Prop) [DecidablePred p] (h_pos : ∀ x ∈ L, p x → 0 ≤ g x) (h_ge : ∀ x ∈ L, p x → f x ≥ g x) : (L.map (fun x => if p x then f x else 0)).sum ≥ (L.map (fun x => if p x then g x else 0)).sum := by
  induction L with
  | nil =>
      simp
  | cons x xs ih =>
      -- restricted hypotheses for the tail
      have h_pos_xs : ∀ y ∈ xs, p y → 0 ≤ g y := by
        intro y hy hpy
        exact h_pos y (List.mem_cons_of_mem x hy) hpy
      have h_ge_xs : ∀ y ∈ xs, p y → f y ≥ g y := by
        intro y hy hpy
        exact h_ge y (List.mem_cons_of_mem x hy) hpy
      -- inequality for the tail sums (in ≤ direction)
      have ih' : (List.map (fun x => if p x then g x else 0) xs).sum ≤
          (List.map (fun x => if p x then f x else 0) xs).sum := by
        have ih0 := ih h_pos_xs h_ge_xs
        -- ih0 : (List.map (fun x => if p x then f x else 0) xs).sum ≥
        --        (List.map (fun x => if p x then g x else 0) xs).sum
        simpa [ge_iff_le] using ih0
      -- inequality for the head terms
      have h_head : (if p x then g x else 0) ≤ (if p x then f x else 0) := by
        by_cases hx : p x
        · have hx_ge : g x ≤ f x := by
            have hg := h_ge x (by simp) hx
            -- hg : f x ≥ g x
            simpa [ge_iff_le] using hg
          simp [hx, hx_ge]
        · simp [hx]
      -- combine head and tail using monotonicity of addition
      have h_all : (if p x then g x else 0) + (List.map (fun x => if p x then g x else 0) xs).sum ≤
          (if p x then f x else 0) + (List.map (fun x => if p x then f x else 0) xs).sum := by
        exact add_le_add h_head ih'
      -- rewrite back to sums over lists
      simpa [List.map, ge_iff_le] using h_all

theorem list_sum_map_ge_of_pointwise_ge {α : Type} (L : List α) (f g : α → ℝ)
  (h_nonneg : ∀ x ∈ L, 0 ≤ g x)
  (h_ge : ∀ x ∈ L, f x ≥ g x) :
  (L.map f).sum ≥ (L.map g).sum := by
  let p : α → Prop := fun _ => True
  have dec_p : DecidablePred p := fun _ => Decidable.isTrue True.intro
  have h_pos_axiom : ∀ x ∈ L, p x → 0 ≤ g x := by
    intros x hx _
    exact h_nonneg x hx
  have h_ge_axiom : ∀ x ∈ L, p x → f x ≥ g x := by
    intros x hx _
    exact h_ge x hx
  have res := sum_map_ite_bound L f g p h_pos_axiom h_ge_axiom
  simp [p] at res
  exact res

theorem seqA_large_part_recurrence (r : ℕ) (a b : Fin r → ℕ) (ha : ∀ i, a i ≥ 1)
(ρ : ℝ) (hρ : ρ > 1)
(M0 : ℕ)
(hM0 : ∀ n : ℕ, n ≥ M0 →
    (∑ i : Fin r, (1 : ℝ) / ((a i * n + b i : ℕ) : ℝ)) ≥ ρ * ((1 : ℝ) / (n : ℝ))) :
  ∀ k : ℕ,
    sum_inv_ge M0 (seqA r a b (k + 1)) ≥ ρ * sum_inv_ge M0 (seqA r a b k) := by
  intro k
  let g : ℕ → ℝ := fun x => if M0 ≤ x then (1 : ℝ) / (x : ℝ) else 0
  have h_eq : ∀ L, sum_inv_ge M0 L = (L.map g).sum := by
    intro L
    rw [sum_inv_ge_eq_sum_ite]
  rw [h_eq, h_eq]
  rw [seqA]
  rw [list_bind_sum_map]
  simp only [List.map_map, Function.comp]
  rw [list_sum_finRange_eq_finset_sum]
  rw [← list_sum_finset_sum_comm]
  rw [← list_sum_map_const_mul]
  apply list_sum_map_ge_of_pointwise_ge
  · intro x hx
    dsimp [g]
    by_cases h : M0 ≤ x
    · rw [if_pos h]
      apply mul_nonneg
      · exact le_of_lt (lt_trans zero_lt_one hρ)
      · exact one_div_nonneg.mpr (Nat.cast_nonneg x)
    · rw [if_neg h]
      simp
  · intro x hx
    dsimp [g]
    by_cases h : M0 ≤ x
    · rw [if_pos h]
      have h_child : ∀ i, M0 ≤ a i * x + b i :=
        fun i => child_ge_of_parent_ge (a i) (b i) x M0 (ha i) h
      have h_lhs : (∑ i : Fin r, g (a i * x + b i)) = ∑ i : Fin r, (1 : ℝ) / ((a i * x + b i) : ℝ) := by
        apply Finset.sum_congr rfl
        intro i _
        dsimp [g]
        rw [if_pos (h_child i)]
        simp only [Nat.cast_add, Nat.cast_mul]
      rw [h_lhs]
      have h_main := hM0 x h
      simp only [Nat.cast_add, Nat.cast_mul] at h_main
      exact h_main
    · rw [if_neg h]
      rw [mul_zero]
      apply Finset.sum_nonneg
      intro i _
      dsimp [g]
      split_ifs
      · exact one_div_nonneg.mpr (Nat.cast_nonneg _)
      · rfl

theorem sum_range_shift_denominators (n : ℕ) : ((List.range n).map (fun x => (1 : ℚ) / ((x + 1 : ℚ)))).sum = ((List.range n).map (fun x => (1 : ℚ) / ((x : ℚ) + 1))).sum := by
  simp

theorem distinct_sum_bound_harmonic (L : List ℕ) (h : L.Nodup) (hpos : ∀ x ∈ L, x > 0) : sum_inv L ≤ harmonic L.length := by
  unfold sum_inv
  calc
    (L.map (fun x => (1 : ℚ) / (x : ℚ))).sum
        ≤ ((List.range L.length).map (fun x => (1 : ℚ) / ((x + 1) : ℚ))).sum :=
          distinct_sum_bound L h hpos
    _ = ((List.range L.length).map (fun x => (1 : ℚ) / ((x : ℚ) + 1))).sum := by
          simpa using (sum_range_shift_denominators L.length)
    _ = harmonic L.length := by
          simpa using (harmonic_list_sum_Q L.length)

theorem seqA_harmonic_bound (r : ℕ) (a b : Fin r → ℕ) (k : ℕ) (h_nodup : (seqA r a b k).Nodup) (h_pos : ∀ x ∈ seqA r a b k, x > 0) : sum_inv (seqA r a b k) ≤ harmonic (r ^ k) := by
  rw [← seqA_length r a b k]
  apply distinct_sum_bound_harmonic
  · exact h_nodup
  · exact h_pos

lemma one_div_le_one_of_one_le (x : ℚ) (hx1 : (1 : ℚ) ≤ x) : 1 / x ≤ (1 : ℚ) := by
  have hx0 : (0 : ℚ) < x := lt_of_lt_of_le (zero_lt_one : (0:ℚ) < 1) hx1
  have hx0' : (0 : ℚ) ≤ x⁻¹ := le_of_lt (inv_pos.mpr hx0)
  have h := mul_le_mul_of_nonneg_right hx1 hx0'
  simpa [one_div, hx0.ne'] using h

theorem sum_reciprocal_le_card {r : ℕ} (a : Fin r → ℕ) (ha : ∀ i, a i ≥ 1) : (∑ i : Fin r, (1 : ℚ) / (a i : ℚ)) ≤ (r : ℚ) := by
  classical
  have hbound : ∀ i : Fin r, (1 : ℚ) / (a i : ℚ) ≤ 1 := by
    intro i
    have hi : (1 : ℚ) ≤ (a i : ℚ) := by
      exact_mod_cast (ha i)
    simpa using (one_div_le_one_of_one_le (a i : ℚ) hi)
  have hs : (∑ i : Fin r, (1 : ℚ) / (a i : ℚ)) ≤ ∑ i : Fin r, (1 : ℚ) := by
    apply Finset.sum_le_sum
    intro i _
    exact hbound i
  have hs' : (∑ i : Fin r, (1 : ℚ) / (a i : ℚ)) ≤ (r : ℚ) := by
    simpa using hs
  exact hs'

lemma sum_filter_inv_le_sum_inv (M : ℕ) (L : List ℕ) :
  ((L.filter (fun x => M ≤ x)).map (fun x => (1 : ℝ) / (x : ℝ))).sum ≤ (L.map (fun x => (1 : ℝ) / (x : ℝ))).sum := by
  induction L with
  | nil => simp
  | cons x xs ih =>
    simp only [List.filter_cons, List.map_cons, List.sum_cons]
    split_ifs
    · exact add_le_add_left ih _
    · apply le_trans ih
      apply le_add_of_nonneg_left
      simp only [one_div, inv_nonneg, Nat.cast_nonneg]

theorem main_theorem (r : ℕ) (a b : Fin r → ℕ) (ha : ∀ i, a i ≥ 1) (hsum : ∑ i : Fin r, (1 : ℚ) / (a i : ℚ) > 1) :
  ∃ k, ¬ (seqA r a b k).Nodup := by
  -- 1. r ≥ 2
  have hr : 2 ≤ r := by
    have h_le : (∑ i : Fin r, (1 : ℚ) / (a i : ℚ)) ≤ r := sum_reciprocal_le_card a ha
    have h_one_lt : 1 < (r : ℚ) := lt_of_lt_of_le hsum h_le
    norm_cast at h_one_lt

  -- 2. Contradiction setup
  by_contra h_nodup
  push_neg at h_nodup

  -- 3. Recurrence parameters
  obtain ⟨ρ, hρ, M0, _, h_rec_point⟩ := per_point_recurrence_large_x r a b ha hsum

  -- 4. Eventually positive
  obtain ⟨K, hK_pos⟩ := seqA_large_part_eventually_positive r a b ha hr h_nodup M0

  -- 5. Define S_k and establish growth
  let S := fun k => sum_inv_ge M0 (seqA r a b k)
  have h_S_growth : ∀ k, S (k + 1) ≥ ρ * S k := by
    intro k
    apply seqA_large_part_recurrence r a b ha ρ hρ M0 h_rec_point

  have h_S_lower : ∀ n, S (K + n) ≥ S K * ρ ^ n := by
    intro n
    induction n with
    | zero => simp
    | succ n ih =>
      rw [pow_succ]
      calc S (K + n + 1)
        _ ≥ ρ * S (K + n) := h_S_growth (K + n)
        _ ≥ ρ * (S K * ρ ^ n) := mul_le_mul_of_nonneg_left ih (le_of_lt (lt_trans zero_lt_one hρ))
        _ = S K * (ρ ^ n * ρ) := by ring

  -- 6. Upper bound
  have h_S_upper : ∀ m, S m ≤ 1 + Real.log ((r : ℝ)^m) := by
    intro m
    let Lm := seqA r a b m
    have h_nodup_m : Lm.Nodup := h_nodup m
    have h_pos_m : ∀ x ∈ Lm, x > 0 := seqA_pos r a b ha m

    have h1 : sum_inv Lm ≤ harmonic (r ^ m) := seqA_harmonic_bound r a b m h_nodup_m h_pos_m
    have h2 : (harmonic (r ^ m) : ℝ) ≤ 1 + Real.log ((r : ℝ)^m) := by
      have := harmonic_Q_le_one_add_log (r^m) (one_le_pow_of_two_le r m hr)
      rwa [Nat.cast_pow] at this

    have h3 : S m ≤ ((sum_inv Lm : ℚ) : ℝ) := sum_inv_ge_le_sum_inv M0 Lm h_pos_m

    calc S m ≤ ((sum_inv Lm : ℚ) : ℝ) := h3
      _ ≤ (harmonic (r^m) : ℝ) := by exact_mod_cast h1
      _ ≤ 1 + Real.log ((r : ℝ)^m) := h2

  -- 7. Contradiction
  let C := S K
  let A := Real.log r
  let B := 1 + (K : ℝ) * A

  have hC_pos : C > 0 := hK_pos

  -- We need to massage the upper bound into A*n + B form
  have h_upper_linear : ∀ n, S (K + n) ≤ A * n + B := by
    intro n
    specialize h_S_upper (K + n)
    rw [log_nat_pow r (K+n) hr] at h_S_upper
    have h_eq : 1 + ((K + n) : ℝ) * Real.log r = A * n + B := by
      dsimp [A, B]
      ring
    push_cast at h_S_upper
    rwa [h_eq] at h_S_upper

  -- exp_dominates_linear C ρ A B -> exists n
  obtain ⟨n, hn⟩ := exp_dominates_linear C hC_pos ρ hρ A B

  have h_lower := h_S_lower n
  have h_upper := h_upper_linear n

  have : C * ρ ^ n ≤ A * n + B := le_trans h_lower h_upper
  linarith