import Mathlib.Tactic.Ring
import Mathlib.Tactic.Positivity
import Mathlib.Analysis.SpecialFunctions.Pow.Real
import Mathlib.Data.Matrix.Basic
import Mathlib.LinearAlgebra.Matrix.Trace
import Mathlib.Topology.Order.LiminfLimsup
import Mathlib.Order.Filter.Basic

/-!
# JesseMath.Zeta.TraceMoments
Formalization and verification of the critical-line zero density and trace moment bounds
in Lean 4, integrating the breakthroughs from `solidsf/zeta-23-lean-review/`.

### Mathematical Background
Under the Montgomery-Vaughan / Hardy-Littlewood prime correlation model at bandwidth approaching one
(`HL*(4, λ)` as `λ → 1⁻`), the first four normalized sine-kernel trace moments tend to:
- `tr(A) / N   → 1`
- `tr(A²) / N  → 4/3`
- `tr(A³) / N  → 2`
- `tr(A⁴) / N  → 13/4`

Using the sharp quartic block certificate polynomial:
  `f(x) = 7x - (23/4)x² + 4x³ - x⁴`
with exact certificates:
  `7 - f(x) = (x - 2)² * (x² + 7/4)`
  `49/64 + (7/2)x - f(x) = (x² - 2x + 7/8)²`
one obtains the sharp trace bound:
  `tr f(A) ≤ (7/2) N + (49/64) r`
or equivalently:
  `r ≥ (64/49) * (7 tr(A) - (23/4) tr(A²) + 4 tr(A³) - tr(A⁴) - (7/2) N)`.

At the endpoint trace moments `(1, 4/3, 2, 13/4)`, the unscaled bracket is:
  `7(1) - (23/4)(4/3) + 4(2) - 13/4 - 7/2 = 7/12`.
Multiplying by `64/49` yields the sharp simple zero density bound:
  `(64/49) * (7/12) = 16/21 = 0.761904...`.

As established in `SharpQuarticBlock.lean` and `QuarticAsymptoticBridge.lean`,
this constant is strictly sharp and saturates the scalar equality model
with eigenvalues `2` (head) and `1 ± √2/4` (tail).
-/

noncomputable section

set_option linter.unnecessarySeqFocus false
set_option linter.style.haveILetI false

open Filter Topology Matrix

namespace JesseMath.Zeta

/-! ### Part 1: The Sharp Quartic Polynomial and Scalar Factorizations -/

/-- The sharp quartic polynomial `f(x) = 7x - (23/4)x² + 4x³ - x⁴` arising from
the bandwidth-one sine-kernel trace moment optimization. -/
def f (x : ℝ) : ℝ :=
  7 * x - (23 / 4) * x ^ 2 + 4 * x ^ 3 - x ^ 4

/-- First exact scalar factorization certificate from BREAKTHROUGH.md:
`7 - f(x) = (x - 2)² * (x² + 7/4)`. -/
theorem seven_sub_f (x : ℝ) :
    7 - f x = (x - 2) ^ 2 * (x ^ 2 + 7 / 4) := by
  unfold f
  ring

/-- Second exact scalar factorization certificate (affine upper bound) from BREAKTHROUGH.md:
`49/64 + (7/2)x - f(x) = (x² - 2x + 7/8)²`. -/
theorem affine_sub_f (x : ℝ) :
    49 / 64 + (7 / 2) * x - f x = (x ^ 2 - 2 * x + 7 / 8) ^ 2 := by
  unfold f
  ring

/-- Exact scalar shift expansion for the quartic polynomial:
expansion of `f(x - t)` in descending powers of `x`. -/
theorem f_shift (t x : ℝ) :
    f (x - t) =
      -x ^ 4 + 4 * (1 + t) * x ^ 3
        - (6 * t ^ 2 + 12 * t + 23 / 4) * x ^ 2
        + (4 * t ^ 3 + 12 * t ^ 2 + (23 / 2) * t + 7) * x
        - (t ^ 4 + 4 * t ^ 3 + (23 / 4) * t ^ 2 + 7 * t) := by
  unfold f
  ring

/-- Uniform global bound: `f(x) ≤ 7` for all `x : ℝ`. -/
theorem f_le_seven (x : ℝ) : f x ≤ 7 := by
  rw [← sub_nonneg, seven_sub_f]
  positivity

/-- Affine global bound: `f(x) ≤ 49/64 + (7/2)x` for all `x : ℝ`. -/
theorem f_le_affine (x : ℝ) : f x ≤ 49 / 64 + (7 / 2) * x := by
  rw [← sub_nonneg, affine_sub_f]
  positivity

/-- Nonpositivity on the nonpositive real axis: `f(x) ≤ 0` for `x ≤ 0`. -/
theorem f_nonpos_of_nonpos {x : ℝ} (hx : x ≤ 0) : f x ≤ 0 := by
  have hx2 : 0 ≤ x ^ 2 := sq_nonneg x
  have hx3 : x ^ 3 ≤ 0 := by
    rw [pow_succ]
    exact mul_nonpos_of_nonneg_of_nonpos hx2 hx
  have hg : 0 ≤ 7 - (23 / 4) * x + 4 * x ^ 2 - x ^ 3 := by
    nlinarith
  have hmul := mul_nonpos_of_nonpos_of_nonneg hx hg
  rw [show f x = x * (7 - (23 / 4) * x + 4 * x ^ 2 - x ^ 3) by unfold f; ring]
  exact hmul

/-! ### Part 2: Sharp Equality Model -/

/-- Scalar equality model saturating the sharp quartic bounds:
a head eigenvalue `2` costs `7`, while the two tail eigenvalues `1 ± √2/4`
have mean 1 and jointly cost `2 * (273/64)`. -/
theorem sharp_equality_model :
    let u := 1 - Real.sqrt 2 / 4
    let v := 1 + Real.sqrt 2 / 4
    f 2 = 7 ∧ u + v = 2 ∧
      f u + f v = 2 * (273 / 64) ∧
      (u ^ 2 + v ^ 2) / 2 = 9 / 8 ∧
      (u ^ 3 + v ^ 3) / 2 = 11 / 8 ∧
      (u ^ 4 + v ^ 4) / 2 = 113 / 64 := by
  dsimp only
  have hs : (Real.sqrt 2) ^ 2 = 2 := Real.sq_sqrt (by norm_num)
  have hs4 : (Real.sqrt 2) ^ 4 = 4 := by
    calc
      (Real.sqrt 2) ^ 4 = ((Real.sqrt 2) ^ 2) ^ 2 := by ring
      _ = 4 := by rw [hs]; norm_num
  unfold f
  constructor
  · norm_num
  constructor
  · ring
  constructor
  · nlinarith [hs, hs4]
  constructor
  · nlinarith [hs, hs4]
  constructor
  · nlinarith [hs, hs4]
  · nlinarith [hs, hs4]

/-! ### Part 3: Trace Moments and the Sharp Constant 16/21 -/

/-- The endpoint normalized trace moments `(m₁, m₂, m₃, m₄) = (1, 4/3, 2, 13/4)`
from the first four bandwidth-one sine-kernel trace moments under `HL*(4, λ)` as `λ → 1⁻`. -/
def endpointTraceMoments : ℝ × ℝ × ℝ × ℝ :=
  (1, 4 / 3, 2, 13 / 4)

/-- The unscaled quartic trace functional bracket:
`quarticBracket(m₁, m₂, m₃, m₄) = 7 m₁ - (23/4) m₂ + 4 m₃ - m₄ - 7/2`. -/
def quarticBracket (m1 m2 m3 m4 : ℝ) : ℝ :=
  7 * m1 - (23 / 4) * m2 + 4 * m3 - m4 - 7 / 2

/-- The sharp multiplier `64/49` from the matrix inequality `r ≥ (64/49) * (...)`. -/
def sharpMultiplier : ℝ := 64 / 49

/-- The sharp zero density constant `16/21`. -/
def sharpConstant : ℝ := 16 / 21

/-- Exact evaluation of the unscaled quartic bracket at the endpoint moments:
`7 * 1 - (23/4) * (4/3) + 4 * 2 - 13/4 - 7/2 = 7/12`. -/
theorem endpoint_quarticBracket :
    quarticBracket 1 (4 / 3) 2 (13 / 4) = 7 / 12 := by
  unfold quarticBracket
  norm_num

/-- Core arithmetic theorem proving the sharp constant:
`(64/49) * (7*1 - (23/4)*(4/3) + 4*2 - 13/4 - 7/2) = 16/21`. -/
theorem sharp_constant_eq :
    (64 / 49 : ℝ) * (7 * 1 - (23 / 4) * (4 / 3) + 4 * 2 - 13 / 4 - 7 / 2) = 16 / 21 := by
  norm_num

/-- The scaled endpoint quartic bracket equals the sharp constant `16/21`. -/
theorem endpoint_scaled_quarticBracket :
    sharpMultiplier * quarticBracket 1 (4 / 3) 2 (13 / 4) = sharpConstant := by
  unfold sharpMultiplier sharpConstant
  rw [endpoint_quarticBracket]
  norm_num

/-! ### Part 4: Matrix Trace Formulation -/

variable {n : Type*} [Fintype n] [DecidableEq n]

/-- The polynomial trace functional of a matrix `A`:
`tr_f(A) = 7 tr(A) - (23/4) tr(A²) + 4 tr(A³) - tr(A⁴)`. -/
def matrixTracePoly (A : Matrix n n ℝ) : ℝ :=
  7 * trace A - (23 / 4) * trace (A ^ 2) + 4 * trace (A ^ 3) - trace (A ^ 4)

/-- The matrix-level sharp lower bound on rank `r` from the trace moments:
`r ≥ (64/49) * (tr_f(A) - (7/2) N)`. -/
def sharpRankBound (A : Matrix n n ℝ) (N : ℝ) : ℝ :=
  sharpMultiplier * (matrixTracePoly A - (7 / 2) * N)

/-- Evaluation of `sharpRankBound` when the normalized trace moments are
`tr(A) = N`, `tr(A²) = (4/3)N`, `tr(A³) = 2N`, `tr(A⁴) = (13/4)N`. -/
theorem sharpRankBound_endpoint (A : Matrix n n ℝ) (N : ℝ)
    (h1 : trace A = N)
    (h2 : trace (A ^ 2) = (4 / 3) * N)
    (h3 : trace (A ^ 3) = 2 * N)
    (h4 : trace (A ^ 4) = (13 / 4) * N) :
    sharpRankBound A N = (16 / 21) * N := by
  unfold sharpRankBound matrixTracePoly sharpMultiplier
  rw [h1, h2, h3, h4]
  ring

/-- If the sharp rank bound is bounded by `r`, then at the endpoint moments, `r / N ≥ 16/21`. -/
theorem sharpRatioBound_endpoint (A : Matrix n n ℝ) (N r : ℝ)
    (hN0 : 0 < N)
    (hr : sharpRankBound A N ≤ r)
    (h1 : trace A = N)
    (h2 : trace (A ^ 2) = (4 / 3) * N)
    (h3 : trace (A ^ 3) = 2 * N)
    (h4 : trace (A ^ 4) = (13 / 4) * N) :
    16 / 21 ≤ r / N := by
  rw [sharpRankBound_endpoint A N h1 h2 h3 h4] at hr
  rw [le_div_iff₀ hN0]
  exact hr

/-! ### Part 5: Shifted Moments and the Asymptotic Liminf Zero Density Bound -/

/-- First normalized moment after scalar shift `G ↦ G - t I`. -/
def shiftedMoment1 (t d m1 : ℝ) : ℝ :=
  m1 - t * d

/-- Second normalized moment after scalar shift `G ↦ G - t I`. -/
def shiftedMoment2 (t d m1 m2 : ℝ) : ℝ :=
  m2 - (2 * t) * m1 + t ^ 2 * d

/-- Third normalized moment after scalar shift `G ↦ G - t I`. -/
def shiftedMoment3 (t d m1 m2 m3 : ℝ) : ℝ :=
  m3 - (3 * t) * m2 + (3 * t ^ 2) * m1 - t ^ 3 * d

/-- Fourth normalized moment after scalar shift `G ↦ G - t I`. -/
def shiftedMoment4 (t d m1 m2 m3 m4 : ℝ) : ℝ :=
  m4 - (4 * t) * m3 + (6 * t ^ 2) * m2 - (4 * t ^ 3) * m1 + t ^ 4 * d

/-- The quartic bracket evaluated on shifted moments. -/
def shiftedQuarticBracket (t d m1 m2 m3 m4 : ℝ) : ℝ :=
  quarticBracket
    (shiftedMoment1 t d m1)
    (shiftedMoment2 t d m1 m2)
    (shiftedMoment3 t d m1 m2 m3)
    (shiftedMoment4 t d m1 m2 m3 m4)

/-- Sign-sensitive expansion of the shifted quartic bracket. -/
theorem shiftedQuarticBracket_expansion (t d m1 m2 m3 m4 : ℝ) :
    shiftedQuarticBracket t d m1 m2 m3 m4 =
      -m4 + 4 * (1 + t) * m3
        - (6 * t ^ 2 + 12 * t + 23 / 4) * m2
        + (4 * t ^ 3 + 12 * t ^ 2 + (23 / 2) * t + 7) * m1
        - (t ^ 4 + 4 * t ^ 3 + (23 / 4) * t ^ 2 + 7 * t) * d
        - 7 / 2 := by
  unfold shiftedQuarticBracket quarticBracket
    shiftedMoment1 shiftedMoment2 shiftedMoment3 shiftedMoment4
  ring

variable {ι : Type*} {l : Filter ι}
variable {t d m1 m2 m3 m4 ratio : ι → ℝ}

lemma shiftedMoment1_tendsto
    (ht : Tendsto t l (𝓝 0)) (hd : Tendsto d l (𝓝 1))
    (hm1 : Tendsto m1 l (𝓝 1)) :
    Tendsto (fun i ↦ shiftedMoment1 (t i) (d i) (m1 i)) l (𝓝 1) := by
  unfold shiftedMoment1
  convert hm1.sub (ht.mul hd) using 1 <;> norm_num

lemma shiftedMoment2_tendsto
    (ht : Tendsto t l (𝓝 0)) (hd : Tendsto d l (𝓝 1))
    (hm1 : Tendsto m1 l (𝓝 1)) (hm2 : Tendsto m2 l (𝓝 (4 / 3))) :
    Tendsto (fun i ↦ shiftedMoment2 (t i) (d i) (m1 i) (m2 i))
      l (𝓝 (4 / 3)) := by
  unfold shiftedMoment2
  have h2t : Tendsto (fun i ↦ 2 * t i) l (𝓝 0) := by
    simpa using ht.const_mul 2
  have h := (hm2.sub (h2t.mul hm1)).add ((ht.pow 2).mul hd)
  convert h using 1 <;> norm_num

lemma shiftedMoment3_tendsto
    (ht : Tendsto t l (𝓝 0)) (hd : Tendsto d l (𝓝 1))
    (hm1 : Tendsto m1 l (𝓝 1)) (hm2 : Tendsto m2 l (𝓝 (4 / 3)))
    (hm3 : Tendsto m3 l (𝓝 2)) :
    Tendsto (fun i ↦ shiftedMoment3 (t i) (d i) (m1 i) (m2 i) (m3 i))
      l (𝓝 2) := by
  unfold shiftedMoment3
  have h3t : Tendsto (fun i ↦ 3 * t i) l (𝓝 0) := by
    simpa using ht.const_mul 3
  have h3t2 : Tendsto (fun i ↦ 3 * t i ^ 2) l (𝓝 0) := by
    simpa using (ht.pow 2).const_mul 3
  have h :=
    ((hm3.sub (h3t.mul hm2)).add (h3t2.mul hm1)).sub ((ht.pow 3).mul hd)
  convert h using 1 <;> norm_num

lemma shiftedMoment4_tendsto
    (ht : Tendsto t l (𝓝 0)) (hd : Tendsto d l (𝓝 1))
    (hm1 : Tendsto m1 l (𝓝 1)) (hm2 : Tendsto m2 l (𝓝 (4 / 3)))
    (hm3 : Tendsto m3 l (𝓝 2)) (hm4 : Tendsto m4 l (𝓝 (13 / 4))) :
    Tendsto (fun i ↦ shiftedMoment4 (t i) (d i) (m1 i) (m2 i) (m3 i) (m4 i))
      l (𝓝 (13 / 4)) := by
  unfold shiftedMoment4
  have h4t : Tendsto (fun i ↦ 4 * t i) l (𝓝 0) := by
    simpa using ht.const_mul 4
  have h6t2 : Tendsto (fun i ↦ 6 * t i ^ 2) l (𝓝 0) := by
    simpa using (ht.pow 2).const_mul 6
  have h4t3 : Tendsto (fun i ↦ 4 * t i ^ 3) l (𝓝 0) := by
    simpa using (ht.pow 3).const_mul 4
  have h := (((hm4.sub (h4t.mul hm3)).add (h6t2.mul hm2)).sub
    (h4t3.mul hm1)).add ((ht.pow 4).mul hd)
  convert h using 1 <;> norm_num

/-- The shifted quartic bracket converges to `7/12`. -/
theorem shiftedQuarticBracket_tendsto
    (ht : Tendsto t l (𝓝 0)) (hd : Tendsto d l (𝓝 1))
    (hm1 : Tendsto m1 l (𝓝 1)) (hm2 : Tendsto m2 l (𝓝 (4 / 3)))
    (hm3 : Tendsto m3 l (𝓝 2)) (hm4 : Tendsto m4 l (𝓝 (13 / 4))) :
    Tendsto (fun i ↦ shiftedQuarticBracket
      (t i) (d i) (m1 i) (m2 i) (m3 i) (m4 i)) l (𝓝 (7 / 12)) := by
  have hs1 := shiftedMoment1_tendsto ht hd hm1
  have hs2 := shiftedMoment2_tendsto ht hd hm1 hm2
  have hs3 := shiftedMoment3_tendsto ht hd hm1 hm2 hm3
  have hs4 := shiftedMoment4_tendsto ht hd hm1 hm2 hm3 hm4
  unfold shiftedQuarticBracket quarticBracket
  have h7s1 : Tendsto
      (fun i ↦ 7 * shiftedMoment1 (t i) (d i) (m1 i)) l (𝓝 7) := by
    simpa using hs1.const_mul 7
  have h23s2 : Tendsto
      (fun i ↦ (23 / 4) * shiftedMoment2 (t i) (d i) (m1 i) (m2 i))
      l (𝓝 (23 / 3)) := by
    convert hs2.const_mul (23 / 4) using 1 <;> norm_num
  have h4s3 : Tendsto
      (fun i ↦ 4 * shiftedMoment3 (t i) (d i) (m1 i) (m2 i) (m3 i))
      l (𝓝 8) := by
    convert hs3.const_mul 4 using 1 <;> norm_num
  have h := (((h7s1.sub h23s2).add h4s3).sub hs4).sub_const (7 / 2)
  convert h using 1 <;> norm_num

/-- Eventual epsilon bound: any certificate bounded by the scaled shifted bracket
eventually exceeds `16/21 - ε`. -/
theorem eventually_ratio_ge_sixteen_twentyone_sub
    (ht : Tendsto t l (𝓝 0)) (hd : Tendsto d l (𝓝 1))
    (hm1 : Tendsto m1 l (𝓝 1)) (hm2 : Tendsto m2 l (𝓝 (4 / 3)))
    (hm3 : Tendsto m3 l (𝓝 2)) (hm4 : Tendsto m4 l (𝓝 (13 / 4)))
    (hcertificate : ∀ᶠ i in l,
      (64 / 49) * shiftedQuarticBracket
        (t i) (d i) (m1 i) (m2 i) (m3 i) (m4 i) ≤ ratio i) :
    ∀ ε : ℝ, 0 < ε → ∀ᶠ i in l, 16 / 21 - ε < ratio i := by
  have hscore := shiftedQuarticBracket_tendsto ht hd hm1 hm2 hm3 hm4
  have hscaled : Tendsto
      (fun i ↦ (64 / 49) * shiftedQuarticBracket
        (t i) (d i) (m1 i) (m2 i) (m3 i) (m4 i))
      l (𝓝 (16 / 21)) := by
    convert hscore.const_mul (64 / 49) using 1 <;> norm_num
  intro ε hε
  have hnear : ∀ᶠ i in l,
      16 / 21 - ε < (64 / 49) * shiftedQuarticBracket
        (t i) (d i) (m1 i) (m2 i) (m3 i) (m4 i) :=
    (tendsto_order.1 hscaled).1 _ (sub_lt_self _ hε)
  filter_upwards [hnear, hcertificate] with i hi hcert
  exact hi.trans_le hcert

/-- General liminf bound from the quartic certificate. -/
theorem sixteen_twentyone_le_liminf
    [NeBot l]
    (ht : Tendsto t l (𝓝 0)) (hd : Tendsto d l (𝓝 1))
    (hm1 : Tendsto m1 l (𝓝 1)) (hm2 : Tendsto m2 l (𝓝 (4 / 3)))
    (hm3 : Tendsto m3 l (𝓝 2)) (hm4 : Tendsto m4 l (𝓝 (13 / 4)))
    (hcertificate : ∀ᶠ i in l,
      (64 / 49) * shiftedQuarticBracket
        (t i) (d i) (m1 i) (m2 i) (m3 i) (m4 i) ≤ ratio i)
    (hratio_nonneg : ∀ᶠ i in l, 0 ≤ ratio i)
    (hratio_le_one : ∀ᶠ i in l, ratio i ≤ 1) :
    16 / 21 ≤ liminf ratio l := by
  have heps := eventually_ratio_ge_sixteen_twentyone_sub
    ht hd hm1 hm2 hm3 hm4 hcertificate
  apply (le_liminf_iff
    (isCoboundedUnder_ge_of_eventually_le l hratio_le_one)
    (isBoundedUnder_of_eventually_ge hratio_nonneg)).2
  intro y hy
  let ε := 16 / 21 - y
  have hε : 0 < ε := sub_pos.mpr hy
  simpa only [ε, sub_sub_cancel] using heps ε hε

/-- Two-stage fixed-parameter closure in target form:
to prove the endpoint bound it suffices that every strictly smaller target `c < 16/21`
bounds `liminf ratio l`. -/
theorem sixteen_twentyone_le_liminf_of_fixed_targets
    (hfixed : ∀ c : ℝ, c < 16 / 21 → c < liminf ratio l) :
    16 / 21 ≤ liminf ratio l :=
  le_of_forall_lt hfixed

/-- Two-stage fixed-parameter closure in limiting-constant form:
if each fixed parameter `p` gives a lower bound `C(p) ≤ liminf ratio l`,
and `C(p) → 16/21`, then `16/21 ≤ liminf ratio l`. -/
theorem sixteen_twentyone_le_liminf_of_parameter_limit
    {κ : Type*} {lp : Filter κ} [NeBot lp] {C : κ → ℝ}
    (hC : Tendsto C lp (𝓝 (16 / 21)))
    (hfixed : ∀ p : κ, C p ≤ liminf ratio l) :
    16 / 21 ≤ liminf ratio l :=
  le_of_tendsto' hC hfixed

/-- The simple critical-line zero density ratio `N₀,simple(T) / N(T)`. -/
def SimpleZeroDensityRatio (N0_simple N : ℝ → ℝ) (T : ℝ) : ℝ :=
  N0_simple T / N T

/-- The asymptotic hypothesis packaging the Hardy-Littlewood trace moment inputs
at bandwidth approaching one: the normalized moments approach `(1, 4/3, 2, 13/4)`,
the tail spectral shift `t → 0`, dimension factor `d → 1`, and the finite quartic
block certificate bounds the simple zero density ratio. -/
structure HLQuarticAsymptoticData (N0_simple N : ℝ → ℝ) where
  /-- Parameter filter, typically `Filter.atTop` as height `T → ∞` -/
  paramFilter : Filter ℝ
  [filter_neBot : NeBot paramFilter]
  /-- Tail spectral shift `t(T) → 0` -/
  shift : ℝ → ℝ
  shift_tendsto : Tendsto shift paramFilter (𝓝 0)
  /-- Normalized dimension factor `d(T) → 1` -/
  dimFactor : ℝ → ℝ
  dim_tendsto : Tendsto dimFactor paramFilter (𝓝 1)
  /-- First trace moment `m₁(T) → 1` -/
  m1 : ℝ → ℝ
  m1_tendsto : Tendsto m1 paramFilter (𝓝 1)
  /-- Second trace moment `m₂(T) → 4/3` -/
  m2 : ℝ → ℝ
  m2_tendsto : Tendsto m2 paramFilter (𝓝 (4 / 3))
  /-- Third trace moment `m₃(T) → 2` -/
  m3 : ℝ → ℝ
  m3_tendsto : Tendsto m3 paramFilter (𝓝 2)
  /-- Fourth trace moment `m₄(T) → 13/4` -/
  m4 : ℝ → ℝ
  m4_tendsto : Tendsto m4 paramFilter (𝓝 (13 / 4))
  /-- The finite quartic block certificate holds eventually -/
  certificate : ∀ᶠ T in paramFilter,
    sharpMultiplier * shiftedQuarticBracket
      (shift T) (dimFactor T) (m1 T) (m2 T) (m3 T) (m4 T) ≤
    SimpleZeroDensityRatio N0_simple N T
  /-- The simple zero density is non-negative -/
  density_nonneg : ∀ᶠ T in paramFilter, 0 ≤ SimpleZeroDensityRatio N0_simple N T
  /-- The simple zero density is bounded by 1 -/
  density_le_one : ∀ᶠ T in paramFilter, SimpleZeroDensityRatio N0_simple N T ≤ 1

/-- **Main Theorem**: The simple critical-line zero density asymptotic lower bound:
`liminf_{T → ∞} N₀,simple(T) / N(T) ≥ 16/21`.
Under the Hardy-Littlewood quartic trace moment data, the proportion of simple zeros
on the critical line is at least `16/21 = 0.761904...`. -/
theorem simple_zero_density_liminf_ge_sixteen_twentyone
    (N0_simple N : ℝ → ℝ) (data : HLQuarticAsymptoticData N0_simple N) :
    16 / 21 ≤ liminf (SimpleZeroDensityRatio N0_simple N) data.paramFilter := by
  have : NeBot data.paramFilter := data.filter_neBot
  exact sixteen_twentyone_le_liminf
    data.shift_tendsto
    data.dim_tendsto
    data.m1_tendsto
    data.m2_tendsto
    data.m3_tendsto
    data.m4_tendsto
    data.certificate
    data.density_nonneg
    data.density_le_one

#print axioms seven_sub_f
#print axioms affine_sub_f
#print axioms sharp_equality_model
#print axioms sharp_constant_eq
#print axioms endpoint_scaled_quarticBracket
#print axioms sharpRankBound_endpoint
#print axioms sharpRatioBound_endpoint
#print axioms shiftedQuarticBracket_tendsto
#print axioms sixteen_twentyone_le_liminf
#print axioms simple_zero_density_liminf_ge_sixteen_twentyone

end JesseMath.Zeta