import MaxLeafFormal.LobsterData
import MaxLeafFormal.Antinorm3

/-!
# Uniform lobsters: the certified inequalities

The two matrix inequalities that `ancillary/lobster_rate_certificate.py` certifies for
`prop:rst-caterpillar`(c) of the paper, for every `L ≥ 1`:

* `fᵀ M_y^{L-1} u_y ≤ (1+y)^{7L} / 1.71927^{7L}` for `y ∈ [0, 3/4]`;
* `(1+y*)^{7L} ≤ C · 1.71928^{7L} · fᵀ M_{y*}^{L-1} u_{y*}` at `y* = 0.3453`, for a constant `C > 0`;

and that `u_y`, `M_y` are the paper's recursion (`Antinorm3.merge`, `Antinorm3.init`, reused).
`prop:rst-caterpillar`(c) itself is not established here.
-/

namespace MaxLeafFormal.Lobster

open Matrix

noncomputable section

/-- `u_y`. -/
def u (y : ℝ) : Fin 3 → ℝ := ![y / 2 * (y + 1/4)^2, y^2 / 4 + y / 8, 1/64]
/-- `M_y`. -/
def M (y : ℝ) : Matrix (Fin 3) (Fin 3) ℝ :=
  !![2 * u y 0, u y 0, u y 0; u y 1 + u y 2, u y 1, u y 1 / 2; 0, u y 2, u y 2 / 2]
/-- `f = (1, 1, 1/2)`. -/
def f : Fin 3 → ℝ := ![1, 1, 1/2]
/-- A vector as a triple. -/
def toT (P : Fin 3 → ℝ) : ℝ × ℝ × ℝ := (P 0, P 1, P 2)

end

theorem u_eq (y : ℝ) : toT (u y) =
    Antinorm3.merge (Antinorm3.merge (Antinorm3.init y 0) (Antinorm3.init y 2)) (Antinorm3.init y 2) := by
  simp only [toT, u, Antinorm3.merge, Antinorm3.init]
  ext <;> simp <;> ring

theorem M_mulVec (y : ℝ) (P : Fin 3 → ℝ) : toT (M y *ᵥ P) = Antinorm3.merge (toT (u y)) (toT P) := by
  simp [toT, M, u, Antinorm3.merge, Matrix.mulVec, dotProduct, Fin.sum_univ_three]
  refine ⟨?_, ?_, ?_⟩ <;> ring

/-! ### The checker (over `ℚ`) -/

section Check

/-- `u_y` over `ℚ`. -/
def uQ (y : ℚ) : ℚ × ℚ × ℚ := (y / 2 * (y + 1/4)^2, y^2 / 4 + y / 8, 1/64)

/-- One cell `(a, b, w₀, w₁, w₂, λ, c)` with `B = 1.71927`:
`0 ≤ a ≤ b`, `w > 0`, `λ > 0`, `wᵀ M_b ≤ λ wᵀ`, `f ≤ c w`, `c (w · u_b) ≤ λ`, `B⁷ λ ≤ (1+a)⁷`. -/
def cellOK (cl : ℚ × ℚ × ℚ × ℚ × ℚ × ℚ × ℚ) : Bool :=
  let a := cl.1; let b := cl.2.1; let w0 := cl.2.2.1; let w1 := cl.2.2.2.1; let w2 := cl.2.2.2.2.1
  let lam := cl.2.2.2.2.2.1; let c := cl.2.2.2.2.2.2
  let u1 := (uQ b).1; let u2 := (uQ b).2.1; let u3 := (uQ b).2.2
  decide (0 ≤ a) && decide (a ≤ b) && decide (0 < w0) && decide (0 < w1) && decide (0 < w2) &&
    decide (0 < lam) &&
    -- columns of wᵀ M_b, with M_b = [[2u1,u1,u1],[u2+u3,u2,u2/2],[0,u3,u3/2]]
    decide (w0 * (2 * u1) + w1 * (u2 + u3) ≤ lam * w0) &&
    decide (w0 * u1 + w1 * u2 + w2 * u3 ≤ lam * w1) &&
    decide (w0 * u1 + w1 * (u2 / 2) + w2 * (u3 / 2) ≤ lam * w2) &&
    decide (1 ≤ c * w0) && decide (1 ≤ c * w1) && decide (1/2 ≤ c * w2) &&
    decide (c * (w0 * u1 + w1 * u2 + w2 * u3) ≤ lam) &&
    decide ((171927/100000 : ℚ)^7 * lam ≤ (1 + a)^7)

/-- The cells chain from `lo` to `hi`: consecutive cells share endpoints. -/
def chainOK : ℚ → List (ℚ × ℚ × ℚ × ℚ × ℚ × ℚ × ℚ) → ℚ → Bool
  | lo, [], hi => decide (lo = hi)
  | lo, cl :: rest, hi => decide (cl.1 = lo) && chainOK cl.2.1 rest hi

/-- The upper certificate at `y* = 0.3453` with `B = 1.71928`: `v > 0`, `λ' > 0`, `M v ≥ λ' v`,
`(1+y*)⁷ ≤ B⁷ λ'`. -/
def upperOK (v : ℚ × ℚ × ℚ) (lam : ℚ) : Bool :=
  let y : ℚ := 3453/10000
  let u1 := (uQ y).1; let u2 := (uQ y).2.1; let u3 := (uQ y).2.2
  decide (0 < v.1) && decide (0 < v.2.1) && decide (0 < v.2.2) && decide (0 < lam) &&
    decide (lam * v.1 ≤ 2 * u1 * v.1 + u1 * v.2.1 + u1 * v.2.2) &&
    decide (lam * v.2.1 ≤ (u2 + u3) * v.1 + u2 * v.2.1 + u2 / 2 * v.2.2) &&
    decide (lam * v.2.2 ≤ u3 * v.2.1 + u3 / 2 * v.2.2) &&
    decide ((1 + y)^7 ≤ (171928/100000 : ℚ)^7 * lam)

end Check

/-! ### The data, checked in the kernel -/

theorem cells0_ok : cells0.all cellOK = true := by decide +kernel
theorem cells1_ok : cells1.all cellOK = true := by decide +kernel
theorem cells2_ok : cells2.all cellOK = true := by decide +kernel
theorem cells3_ok : cells3.all cellOK = true := by decide +kernel
theorem cells4_ok : cells4.all cellOK = true := by decide +kernel
theorem cells5_ok : cells5.all cellOK = true := by decide +kernel
theorem cells6_ok : cells6.all cellOK = true := by decide +kernel
theorem cells7_ok : cells7.all cellOK = true := by decide +kernel
theorem cells8_ok : cells8.all cellOK = true := by decide +kernel
theorem cells9_ok : cells9.all cellOK = true := by decide +kernel
theorem chain_ok : chainOK 0 chunks.flatten (3/4) = true := by decide +kernel
theorem upper_ok : upperOK vUp lamUp = true := by decide +kernel

/-! ### Soundness and the theorems -/

/-- Coordinates of `A *ᵥ x` for `3 × 3` matrices. -/
theorem mulVec_apply3 (A : Matrix (Fin 3) (Fin 3) ℝ) (x : Fin 3 → ℝ) (i : Fin 3) :
    (A *ᵥ x) i = A i 0 * x 0 + A i 1 * x 1 + A i 2 * x 2 := by
  simp [Matrix.mulVec, dotProduct, Fin.sum_univ_three]

theorem f_dot (x : Fin 3 → ℝ) : f ⬝ᵥ x = x 0 + x 1 + x 2 / 2 := by
  simp [f, dotProduct, Fin.sum_univ_three] <;> ring

theorem M_mulVec_coord (y : ℝ) (X : Fin 3 → ℝ) :
    (M y *ᵥ X) 0 = 2 * u y 0 * X 0 + u y 0 * X 1 + u y 0 * X 2 ∧
    (M y *ᵥ X) 1 = (u y 1 + u y 2) * X 0 + u y 1 * X 1 + u y 1 / 2 * X 2 ∧
    (M y *ᵥ X) 2 = u y 2 * X 1 + u y 2 / 2 * X 2 := by
  refine ⟨?_, ?_, ?_⟩ <;> rw [mulVec_apply3] <;> simp [M] <;> ring

theorem u_nonneg {y : ℝ} (hy : 0 ≤ y) (i : Fin 3) : 0 ≤ u y i := by
  fin_cases i <;> simp [u] <;> nlinarith [mul_nonneg hy (sq_nonneg (y + 1/4)), sq_nonneg y]

theorem u_mono {y b : ℝ} (hy : 0 ≤ y) (hyb : y ≤ b) (i : Fin 3) : u y i ≤ u b i := by
  have hb : 0 ≤ b := le_trans hy hyb
  fin_cases i <;> simp [u] <;>
    nlinarith [mul_nonneg (sub_nonneg.2 hyb) (sq_nonneg (b + 1/4)),
      mul_nonneg (mul_nonneg hy (sub_nonneg.2 hyb)) (show (0:ℝ) ≤ b + y + 1/2 by linarith),
      mul_nonneg (sub_nonneg.2 hyb) (show (0:ℝ) ≤ b + y by linarith)]

theorem M_le_M {y b : ℝ} (hy : 0 ≤ y) (hyb : y ≤ b) (i j : Fin 3) :
    0 ≤ M y i j ∧ M y i j ≤ M b i j := by
  have h0 := u_nonneg hy 0
  have h1 := u_nonneg hy 1
  have h2 := u_nonneg hy 2
  have m0 := u_mono hy hyb 0
  have m1 := u_mono hy hyb 1
  have m2 := u_mono hy hyb 2
  fin_cases i <;> fin_cases j <;> simp [M] <;> (try constructor) <;> linarith

theorem M_nonneg {y : ℝ} (hy : 0 ≤ y) (i j : Fin 3) : 0 ≤ M y i j :=
  (M_le_M hy le_rfl i j).1

theorem mulVec_mono_mat {A B : Matrix (Fin 3) (Fin 3) ℝ} {x x' : Fin 3 → ℝ}
    (_hA : ∀ i j, 0 ≤ A i j) (hAB : ∀ i j, A i j ≤ B i j) (hx : ∀ j, 0 ≤ x j)
    (hxx : ∀ j, x j ≤ x' j) (hB : ∀ i j, 0 ≤ B i j) (i : Fin 3) :
    (A *ᵥ x) i ≤ (B *ᵥ x') i := by
  rw [mulVec_apply3, mulVec_apply3]
  have h0 := mul_le_mul (hAB i 0) (hxx 0) (hx 0) (hB i 0)
  have h1 := mul_le_mul (hAB i 1) (hxx 1) (hx 1) (hB i 1)
  have h2 := mul_le_mul (hAB i 2) (hxx 2) (hx 2) (hB i 2)
  linarith

theorem mulVec_nonneg3 {A : Matrix (Fin 3) (Fin 3) ℝ} {x : Fin 3 → ℝ} (hA : ∀ i j, 0 ≤ A i j)
    (hx : ∀ j, 0 ≤ x j) (i : Fin 3) : 0 ≤ (A *ᵥ x) i := by
  rw [mulVec_apply3]
  have h0 := mul_nonneg (hA i 0) (hx 0)
  have h1 := mul_nonneg (hA i 1) (hx 1)
  have h2 := mul_nonneg (hA i 2) (hx 2)
  linarith

theorem iter_nonneg {y : ℝ} (hy : 0 ≤ y) (n : ℕ) (i : Fin 3) : 0 ≤ (M y ^ n *ᵥ u y) i := by
  induction n generalizing i with
  | zero => simpa using u_nonneg hy i
  | succ n ih =>
    rw [pow_succ' (M y) n, ← mulVec_mulVec]
    exact mulVec_nonneg3 (fun i j => M_nonneg hy i j) (fun j => ih j) i

theorem iter_mono {y b : ℝ} (hy : 0 ≤ y) (hyb : y ≤ b) (n : ℕ) (i : Fin 3) :
    (M y ^ n *ᵥ u y) i ≤ (M b ^ n *ᵥ u b) i := by
  induction n generalizing i with
  | zero => simpa using u_mono hy hyb i
  | succ n ih =>
    rw [pow_succ' (M y) n, pow_succ' (M b) n, ← mulVec_mulVec, ← mulVec_mulVec]
    exact mulVec_mono_mat (fun i j => M_nonneg hy i j) (fun i j => (M_le_M hy hyb i j).2)
      (fun j => iter_nonneg hy n j) ih (fun i j => M_nonneg (le_trans hy hyb) i j) i

/-- The lower bound for one cell, in `ℝ`. -/
theorem lower_core (b lam c w0 w1 w2 : ℝ) (hb : 0 ≤ b) (hlam : 0 ≤ lam) (hc : 0 ≤ c)
    (hcol0 : w0 * (2 * u b 0) + w1 * (u b 1 + u b 2) ≤ lam * w0)
    (hcol1 : w0 * u b 0 + w1 * u b 1 + w2 * u b 2 ≤ lam * w1)
    (hcol2 : w0 * u b 0 + w1 * (u b 1 / 2) + w2 * (u b 2 / 2) ≤ lam * w2)
    (hf0 : 1 ≤ c * w0) (hf1 : 1 ≤ c * w1) (hf2 : 1 / 2 ≤ c * w2)
    (hcu : c * (w0 * u b 0 + w1 * u b 1 + w2 * u b 2) ≤ lam)
    (y : ℝ) (hy : 0 ≤ y) (hyb : y ≤ b) (k : ℕ) :
    f ⬝ᵥ (M y ^ k *ᵥ u y) ≤ lam ^ (k + 1) := by
  have hW : ∀ n : ℕ, w0 * (M b ^ n *ᵥ u b) 0 + w1 * (M b ^ n *ᵥ u b) 1 + w2 * (M b ^ n *ᵥ u b) 2
      ≤ lam ^ n * (w0 * u b 0 + w1 * u b 1 + w2 * u b 2) := by
    intro n
    induction n with
    | zero => simp
    | succ n ih =>
      have X0 := iter_nonneg hb n 0
      have X1 := iter_nonneg hb n 1
      have X2 := iter_nonneg hb n 2
      rw [pow_succ' (M b) n, ← mulVec_mulVec]
      obtain ⟨e0, e1, e2⟩ := M_mulVec_coord b (M b ^ n *ᵥ u b)
      rw [e0, e1, e2, pow_succ lam n]
      have s0 := mul_le_mul_of_nonneg_right hcol0 X0
      have s1 := mul_le_mul_of_nonneg_right hcol1 X1
      have s2 := mul_le_mul_of_nonneg_right hcol2 X2
      have h2 := mul_le_mul_of_nonneg_left ih hlam
      nlinarith
  rw [f_dot, pow_succ lam k]
  have hXb := fun i => iter_nonneg hb k i
  have g1 := mul_le_mul_of_nonneg_right hf0 (hXb 0)
  have g2 := mul_le_mul_of_nonneg_right hf1 (hXb 1)
  have g3 := mul_le_mul_of_nonneg_right hf2 (hXb 2)
  have g4 := mul_le_mul_of_nonneg_left (hW k) hc
  have g5 := mul_le_mul_of_nonneg_left hcu (pow_nonneg hlam k)
  have hm0 := iter_mono hy hyb k 0
  have hm1 := iter_mono hy hyb k 1
  have hm2 := iter_mono hy hyb k 2
  linarith

/-- The upper bound certificate at `y*`, for a vector `x` with `M x ≥ λ x`. -/
theorem upper_core (y lam : ℝ) (x : Fin 3 → ℝ) (hy : 0 ≤ y) (hlam : 0 < lam) (hx : ∀ i, 0 < x i)
    (hu : ∀ i, 0 < u y i) (hM0 : lam * x 0 ≤ (M y *ᵥ x) 0) (hM1 : lam * x 1 ≤ (M y *ᵥ x) 1)
    (hM2 : lam * x 2 ≤ (M y *ᵥ x) 2) :
    ∃ c : ℝ, 0 < c ∧ ∀ k : ℕ, c * lam ^ k * (f ⬝ᵥ x) ≤ f ⬝ᵥ (M y ^ k *ᵥ u y) := by
  obtain ⟨c, hcpos, hcle⟩ : ∃ c : ℝ, 0 < c ∧ ∀ i, c * x i ≤ u y i := by
    refine ⟨min (u y 0 / x 0) (min (u y 1 / x 1) (u y 2 / x 2)),
      lt_min (div_pos (hu 0) (hx 0)) (lt_min (div_pos (hu 1) (hx 1)) (div_pos (hu 2) (hx 2))), ?_⟩
    intro i
    fin_cases i
    · exact (le_div_iff₀ (hx 0)).1 (min_le_left _ _)
    · exact (le_div_iff₀ (hx 1)).1 (le_trans (min_le_right _ _) (min_le_left _ _))
    · exact (le_div_iff₀ (hx 2)).1 (le_trans (min_le_right _ _) (min_le_right _ _))
  have hMv : ∀ i, lam * x i ≤ (M y *ᵥ x) i := by
    intro i
    fin_cases i
    · exact hM0
    · exact hM1
    · exact hM2
  have hind : ∀ k : ℕ, ∀ i, c * lam ^ k * x i ≤ (M y ^ k *ᵥ u y) i := by
    intro k
    induction k with
    | zero => intro i; simpa using hcle i
    | succ n ih =>
      intro i
      have hpos : 0 ≤ c * lam ^ n := mul_nonneg hcpos.le (pow_nonneg hlam.le n)
      rw [pow_succ' (M y) n, ← mulVec_mulVec, pow_succ lam n]
      have h1 := mulVec_mono_mat (A := M y) (B := M y) (x := fun j => c * lam ^ n * x j)
        (x' := M y ^ n *ᵥ u y) (fun i j => M_nonneg hy i j) (fun i j => le_rfl)
        (fun j => mul_nonneg hpos (hx j).le) (fun j => ih j) (fun i j => M_nonneg hy i j) i
      have h2 : (M y *ᵥ (fun j => c * lam ^ n * x j)) i = c * lam ^ n * (M y *ᵥ x) i := by
        simp only [mulVec_apply3]; ring
      have h3 := mul_le_mul_of_nonneg_left (hMv i) hpos
      linarith
  refine ⟨c, hcpos, fun k => ?_⟩
  rw [f_dot, f_dot]
  have a0 := hind k 0
  have a1 := hind k 1
  have a2 := hind k 2
  linarith

theorem uQ_cast (b : ℚ) : ((uQ b).1 : ℝ) = u (b:ℝ) 0 ∧ ((uQ b).2.1 : ℝ) = u (b:ℝ) 1 ∧
    ((uQ b).2.2 : ℝ) = u (b:ℝ) 2 := by
  simp only [uQ, u]
  refine ⟨?_, ?_, ?_⟩ <;> push_cast <;> ring

theorem all_cells_ok : ∀ cl ∈ chunks.flatten, cellOK cl = true := by
  intro cl hcl
  obtain ⟨l, hl, hcl'⟩ := List.mem_flatten.1 hcl
  simp [chunks] at hl
  rcases hl with rfl | rfl | rfl | rfl | rfl | rfl | rfl | rfl | rfl | rfl
  · exact List.all_eq_true.1 cells0_ok cl hcl'
  · exact List.all_eq_true.1 cells1_ok cl hcl'
  · exact List.all_eq_true.1 cells2_ok cl hcl'
  · exact List.all_eq_true.1 cells3_ok cl hcl'
  · exact List.all_eq_true.1 cells4_ok cl hcl'
  · exact List.all_eq_true.1 cells5_ok cl hcl'
  · exact List.all_eq_true.1 cells6_ok cl hcl'
  · exact List.all_eq_true.1 cells7_ok cl hcl'
  · exact List.all_eq_true.1 cells8_ok cl hcl'
  · exact List.all_eq_true.1 cells9_ok cl hcl'

theorem chain_cover (cells : List (ℚ × ℚ × ℚ × ℚ × ℚ × ℚ × ℚ)) : ∀ (lo hi : ℚ), cells ≠ [] →
    chainOK lo cells hi = true → ∀ y : ℝ, (lo:ℝ) ≤ y → y ≤ (hi:ℝ) →
    ∃ cl ∈ cells, (cl.1:ℝ) ≤ y ∧ y ≤ (cl.2.1:ℝ) := by
  induction cells with
  | nil => intro lo hi hne; exact absurd rfl hne
  | cons cl rest ih =>
    intro lo hi _ hch y hlo hhi
    simp only [chainOK, Bool.and_eq_true, decide_eq_true_eq] at hch
    obtain ⟨hl, hrest⟩ := hch
    by_cases hy : y ≤ (cl.2.1 : ℝ)
    · exact ⟨cl, by simp, by rw [hl]; exact hlo, hy⟩
    · push_neg at hy
      by_cases hr : rest = []
      · subst hr
        simp only [chainOK, decide_eq_true_eq] at hrest
        exfalso
        have h1 : (cl.2.1 : ℝ) = (hi : ℝ) := by rw [hrest]
        linarith
      · obtain ⟨cl', hmem, h1, h2⟩ := ih cl.2.1 hi hr hrest y hy.le hhi
        exact ⟨cl', by simp [hmem], h1, h2⟩

theorem upper_real (v : ℚ × ℚ × ℚ) (lam : ℚ) (h : upperOK v lam = true) :
    0 < (v.1:ℝ) ∧ 0 < (v.2.1:ℝ) ∧ 0 < (v.2.2:ℝ) ∧ 0 < (lam:ℝ) ∧
    (lam:ℝ) * v.1 ≤ 2 * u (3453/10000:ℝ) 0 * v.1 + u (3453/10000:ℝ) 0 * v.2.1 + u (3453/10000:ℝ) 0 * v.2.2 ∧
    (lam:ℝ) * v.2.1 ≤ (u (3453/10000:ℝ) 1 + u (3453/10000:ℝ) 2) * v.1 + u (3453/10000:ℝ) 1 * v.2.1
      + u (3453/10000:ℝ) 1 / 2 * v.2.2 ∧
    (lam:ℝ) * v.2.2 ≤ u (3453/10000:ℝ) 2 * v.2.1 + u (3453/10000:ℝ) 2 / 2 * v.2.2 ∧
    (1 + 3453/10000 : ℝ)^7 ≤ (171928/100000 : ℝ)^7 * lam := by
  simp only [upperOK, Bool.and_eq_true, decide_eq_true_eq, and_assoc] at h
  obtain ⟨h1, h2, h3, h4, r1, r2, r3, hB⟩ := h
  have e := uQ_cast (3453/10000 : ℚ)
  push_cast at e
  obtain ⟨e0, e1, e2⟩ := e
  have hB' := (Rat.cast_le (K := ℝ)).2 hB
  push_cast at hB'
  refine ⟨by exact_mod_cast h1, by exact_mod_cast h2, by exact_mod_cast h3, by exact_mod_cast h4,
    ?_, ?_, ?_, by linarith⟩
  · rw [← e0]; exact_mod_cast r1
  · rw [← e1, ← e2]; exact_mod_cast r2
  · rw [← e2]; exact_mod_cast r3

theorem cell_bound (a b w0 w1 w2 lam c : ℚ) (hc : cellOK (a, b, w0, w1, w2, lam, c) = true)
    (L : ℕ) (hL : 1 ≤ L) (y : ℝ) (ha : (a:ℝ) ≤ y) (hby : y ≤ (b:ℝ)) :
    f ⬝ᵥ (M y ^ (L - 1) *ᵥ u y) ≤ (1 + y) ^ (7 * L) / (171927/100000 : ℝ) ^ (7 * L) := by
  simp only [cellOK, Bool.and_eq_true, decide_eq_true_eq, and_assoc] at hc
  obtain ⟨h1, h2, h3, h4, h5, h6, h7, h8, h9, h10, h11, h12, h13, h14⟩ := hc
  have hu := uQ_cast b
  have ha0 : (0:ℝ) ≤ a := by exact_mod_cast h1
  have hb0 : (0:ℝ) ≤ b := le_trans ha0 (by exact_mod_cast h2)
  have hw0 : (0:ℝ) < w0 := by exact_mod_cast h3
  have hlam : (0:ℝ) ≤ lam := by exact_mod_cast h6.le
  have hf0 : (1:ℝ) ≤ (c:ℝ) * (w0:ℝ) := by exact_mod_cast h10
  have hf1 : (1:ℝ) ≤ (c:ℝ) * (w1:ℝ) := by exact_mod_cast h11
  have hf2 : (1/2:ℝ) ≤ (c:ℝ) * (w2:ℝ) := by
    have h := (Rat.cast_le (K := ℝ)).2 h12
    push_cast at h
    linarith
  have hc0 : (0:ℝ) ≤ c := by
    by_contra hneg
    push_neg at hneg
    nlinarith
  have hcol0 : (w0:ℝ) * (2 * u (b:ℝ) 0) + (w1:ℝ) * (u (b:ℝ) 1 + u (b:ℝ) 2) ≤ (lam:ℝ) * w0 := by
    rw [← hu.1, ← hu.2.1, ← hu.2.2]; exact_mod_cast h7
  have hcol1 : (w0:ℝ) * u (b:ℝ) 0 + (w1:ℝ) * u (b:ℝ) 1 + (w2:ℝ) * u (b:ℝ) 2 ≤ (lam:ℝ) * w1 := by
    rw [← hu.1, ← hu.2.1, ← hu.2.2]; exact_mod_cast h8
  have hcol2 : (w0:ℝ) * u (b:ℝ) 0 + (w1:ℝ) * (u (b:ℝ) 1 / 2) + (w2:ℝ) * (u (b:ℝ) 2 / 2)
      ≤ (lam:ℝ) * w2 := by
    rw [← hu.1, ← hu.2.1, ← hu.2.2]; exact_mod_cast h9
  have hcu : (c:ℝ) * ((w0:ℝ) * u (b:ℝ) 0 + (w1:ℝ) * u (b:ℝ) 1 + (w2:ℝ) * u (b:ℝ) 2) ≤ (lam:ℝ) := by
    rw [← hu.1, ← hu.2.1, ← hu.2.2]; exact_mod_cast h13
  have hBr : (171927/100000 : ℝ) ^ 7 * (lam:ℝ) ≤ (1 + (a:ℝ)) ^ 7 := by
    have h := (Rat.cast_le (K := ℝ)).2 h14
    push_cast at h
    linarith
  have key := lower_core (b:ℝ) (lam:ℝ) (c:ℝ) (w0:ℝ) (w1:ℝ) (w2:ℝ) hb0 hlam hc0 hcol0 hcol1
    hcol2 hf0 hf1 hf2 hcu y (le_trans ha0 ha) hby (L - 1)
  rw [Nat.sub_add_cancel hL] at key
  have hlamB : (lam:ℝ) ≤ (1 + y) ^ 7 / (171927/100000 : ℝ) ^ 7 := by
    rw [le_div_iff₀ (by positivity)]
    have h1' : (1 + (a:ℝ)) ^ 7 ≤ (1 + y) ^ 7 := pow_le_pow_left₀ (by linarith) (by linarith) 7
    linarith
  have hdiv : (1 + y) ^ (7 * L) / (171927/100000 : ℝ) ^ (7 * L)
      = ((1 + y) ^ 7 / (171927/100000 : ℝ) ^ 7) ^ L := by
    rw [pow_mul (1 + y) 7 L, pow_mul (171927/100000 : ℝ) 7 L]
    exact (div_pow _ _ _).symm
  rw [hdiv]
  exact le_trans key (pow_le_pow_left₀ hlam hlamB L)

theorem lobster_lower (L : ℕ) (hL : 1 ≤ L) (y : ℝ) (hy : y ∈ Set.Icc (0:ℝ) (3/4)) :
    f ⬝ᵥ (M y ^ (L - 1) *ᵥ u y) ≤ (1 + y) ^ (7 * L) / (171927/100000 : ℝ) ^ (7 * L) := by
  obtain ⟨hy0, hy1⟩ := hy
  have hne : chunks.flatten ≠ [] := by
    intro h
    have hc := chain_ok
    rw [h] at hc
    norm_num [chainOK] at hc
  obtain ⟨⟨a, b, w0, w1, w2, lam, c⟩, hcl, ha, hb⟩ :=
    chain_cover _ 0 (3/4) hne chain_ok y (by push_cast; exact hy0) (by push_cast; exact hy1)
  exact cell_bound a b w0 w1 w2 lam c (all_cells_ok _ hcl) L hL y ha hb

theorem lobster_upper : ∃ C : ℝ, 0 < C ∧ ∀ L : ℕ, 1 ≤ L →
    (1 + 3453/10000 : ℝ) ^ (7 * L)
      ≤ C * (171928/100000 : ℝ) ^ (7 * L) * (f ⬝ᵥ (M (3453/10000) ^ (L - 1) *ᵥ u (3453/10000))) := by
  obtain ⟨hv1, hv2, hv3, hlamR, r1, r2, r3, hB⟩ := upper_real vUp lamUp upper_ok
  have hx : ∀ i, 0 < (![((vUp.1 : ℚ) : ℝ), ((vUp.2.1 : ℚ) : ℝ), ((vUp.2.2 : ℚ) : ℝ)] : Fin 3 → ℝ) i := by
    intro i
    fin_cases i <;> simp [hv1, hv2, hv3]
  have hu : ∀ i, 0 < u (3453/10000 : ℝ) i := by
    intro i
    fin_cases i <;> norm_num [u]
  have hM0 : (lamUp:ℝ) * (![((vUp.1 : ℚ) : ℝ), ((vUp.2.1 : ℚ) : ℝ), ((vUp.2.2 : ℚ) : ℝ)] : Fin 3 → ℝ) 0
      ≤ (M (3453/10000 : ℝ) *ᵥ (![((vUp.1 : ℚ) : ℝ), ((vUp.2.1 : ℚ) : ℝ), ((vUp.2.2 : ℚ) : ℝ)] : Fin 3 → ℝ)) 0 := by
    rw [(M_mulVec_coord _ _).1]; simpa using r1
  have hM1 : (lamUp:ℝ) * (![((vUp.1 : ℚ) : ℝ), ((vUp.2.1 : ℚ) : ℝ), ((vUp.2.2 : ℚ) : ℝ)] : Fin 3 → ℝ) 1
      ≤ (M (3453/10000 : ℝ) *ᵥ (![((vUp.1 : ℚ) : ℝ), ((vUp.2.1 : ℚ) : ℝ), ((vUp.2.2 : ℚ) : ℝ)] : Fin 3 → ℝ)) 1 := by
    rw [(M_mulVec_coord _ _).2.1]; simpa using r2
  have hM2 : (lamUp:ℝ) * (![((vUp.1 : ℚ) : ℝ), ((vUp.2.1 : ℚ) : ℝ), ((vUp.2.2 : ℚ) : ℝ)] : Fin 3 → ℝ) 2
      ≤ (M (3453/10000 : ℝ) *ᵥ (![((vUp.1 : ℚ) : ℝ), ((vUp.2.1 : ℚ) : ℝ), ((vUp.2.2 : ℚ) : ℝ)] : Fin 3 → ℝ)) 2 := by
    rw [(M_mulVec_coord _ _).2.2]; simpa using r3
  obtain ⟨c, hc, hk⟩ := upper_core (3453/10000 : ℝ) (lamUp:ℝ)
    (![((vUp.1 : ℚ) : ℝ), ((vUp.2.1 : ℚ) : ℝ), ((vUp.2.2 : ℚ) : ℝ)] : Fin 3 → ℝ)
    (by norm_num) hlamR hx hu hM0 hM1 hM2
  obtain ⟨S, hSdef⟩ : ∃ S : ℝ, S = ((vUp.1 : ℚ) : ℝ) + ((vUp.2.1 : ℚ) : ℝ) + ((vUp.2.2 : ℚ) : ℝ) / 2 :=
    ⟨_, rfl⟩
  have hSv : f ⬝ᵥ (![((vUp.1 : ℚ) : ℝ), ((vUp.2.1 : ℚ) : ℝ), ((vUp.2.2 : ℚ) : ℝ)] : Fin 3 → ℝ) = S := by
    rw [f_dot, hSdef]; simp
  rw [hSv] at hk
  have hSpos : 0 < S := by
    rw [hSdef]; linarith
  refine ⟨(lamUp:ℝ) / (c * S), div_pos hlamR (mul_pos hc hSpos), fun L hL => ?_⟩
  obtain ⟨k, rfl⟩ : ∃ k, L = k + 1 := ⟨L - 1, by omega⟩
  rw [Nat.add_sub_cancel]
  have hCS : (lamUp:ℝ) / (c * S) * (c * (lamUp:ℝ) ^ k * S) = (lamUp:ℝ) ^ (k + 1) := by
    rw [div_mul_eq_mul_div, div_eq_iff (mul_pos hc hSpos).ne']
    ring
  have hLk : (lamUp:ℝ) ^ (k + 1) ≤ (lamUp:ℝ) / (c * S) *
      (f ⬝ᵥ (M (3453/10000 : ℝ) ^ k *ᵥ u (3453/10000 : ℝ))) := by
    rw [← hCS]
    exact mul_le_mul_of_nonneg_left (hk k) (div_pos hlamR (mul_pos hc hSpos)).le
  have hpow : (1 + 3453/10000 : ℝ) ^ (7 * (k + 1))
      ≤ (171928/100000 : ℝ) ^ (7 * (k + 1)) * (lamUp:ℝ) ^ (k + 1) := by
    calc (1 + 3453/10000 : ℝ) ^ (7 * (k + 1)) = ((1 + 3453/10000 : ℝ) ^ 7) ^ (k + 1) := pow_mul _ _ _
      _ ≤ ((171928/100000 : ℝ) ^ 7 * (lamUp:ℝ)) ^ (k + 1) :=
          pow_le_pow_left₀ (by positivity) hB (k + 1)
      _ = (171928/100000 : ℝ) ^ (7 * (k + 1)) * (lamUp:ℝ) ^ (k + 1) := by
          rw [mul_pow, ← pow_mul]
  calc _ ≤ _ := hpow
    _ ≤ (171928/100000 : ℝ) ^ (7 * (k + 1)) * ((lamUp:ℝ) / (c * S) *
          (f ⬝ᵥ (M (3453/10000 : ℝ) ^ k *ᵥ u (3453/10000 : ℝ)))) :=
          mul_le_mul_of_nonneg_left hLk (by positivity)
    _ = _ := by ring

end MaxLeafFormal.Lobster