Round-six κ development, definitions part 60

-- Round-six κ development (integer-multiplication-kappa), generated by skeleton subtraction.

import Definitions.Def_IntMul_Kappa_K_Defs33
import Definitions.Def_IntMul_Kappa_R6_Defs1
import Mathlib
import Definitions.Def_IntMul_HvdH_M_Defs8
import Definitions.Def_IntMul_MultitapeModel
import Definitions.Def_IntMul_HvdH_Numeric
import Definitions.Def_IntMul_HvdH_ExpApprox
import Definitions.Def_IntMul_HvdH_Twiddle
import Definitions.Def_IntMul_HvdH_RealExp
import Definitions.Def_IntMul_HvdH_SApprox
import Definitions.Def_IntMul_HvdH_EApprox
import Definitions.Def_IntMul_HvdH_JApprox
import Definitions.Def_IntMul_HvdH_M_Defs1
import Definitions.Def_IntMul_Kappa_R6_Defs10
import Definitions.Def_IntMul_Kappa_R6_Defs11
import Definitions.Def_IntMul_HvdH_M_Defs2
import Definitions.Def_IntMul_HvdH_M_Defs5
import Definitions.Def_IntMul_Kappa_R6_Defs3
import Definitions.Def_IntMul_Kappa_R6_Defs23
import Definitions.Def_IntMul_HvdH_M_Defs3
import Definitions.Def_IntMul_Kappa_R6_Defs16
import Definitions.Def_IntMul_Kappa_K_Defs10
import Definitions.Def_IntMul_Kappa_K_Defs11
import Definitions.Def_IntMul_Kappa_K_Defs12
import Definitions.Def_IntMul_Kappa_R6_Defs12
import Definitions.Def_IntMul_Kappa_K_Defs13
import Definitions.Def_IntMul_Kappa_R6_Defs14
import Definitions.Def_IntMul_Kappa_K_Defs16
import Definitions.Def_IntMul_Kappa_R6_Defs13
import Definitions.Def_IntMul_Kappa_K_Defs9
import Definitions.Def_IntMul_Kappa_K_Defs23
import Definitions.Def_IntMul_Kappa_K_Defs21

set_option autoImplicit false
set_option relaxedAutoImplicit false
set_option maxHeartbeats 1000000

section
-- Source module: NLogN.Flat
/-!
# Flat code with jumps, and compilation of structured programs into it

Instructions carry absolute addresses. `ret` with an empty control stack halts.
`exec_flat`: a structured run of cost `c` is `c` flat steps.
-/

namespace NLogN

variable {S P : ℕ}

/-! ## Layout -/

theorem Stmt.size_pos (s : Stmt S P) : 0 < s.size := by
  induction s <;> simp [Stmt.size] <;> omega

/-! ## A concrete layout -/

end NLogN

end

section
-- Source module: NLogN.SSBase
/-!
# Schönhage–Strassen in the item language: stack layout and Fermat-ring arithmetic

The program uses `NS` item stacks with fixed roles, so that stack (in)equalities are decided by
`simp`.  The current ring `ℤ/(2^n+1)` is given by the ruler `2^n` on top of `RN` and the modulus
`2^n+1` on top of `RF`.
-/

namespace NLogN

open IntMul

/-! ## Fermat reduction -/

variable {P : ℕ}

end NLogN

namespace NLogN

open IntMul

variable {P : ℕ} {I : List Sym} {iprog : Fin P → IStmt NS P}

theorem zsz_le_of_abs {z : ℤ} {m : ℕ} (h : |z| < 2 ^ m) : zsz z ≤ m := by
  apply zsz_le_of_lt
  have : ((z.natAbs : ℕ) : ℤ) < ((2 ^ m : ℕ) : ℤ) := by push_cast; exact h
  exact_mod_cast this

end NLogN

end

section
-- Source module: NLogN.HReg
/-!
# Register programs

`RS` is a small structured language over the register stacks `0 … 55`: assignments of
expressions, pushes and pops, conditionals and counted loops.  `RS.run` is its functional
semantics on `IState`s, and `RS.code_runs` shows that the compiled `IStmt` realises it with cost
`RS.cost`.
-/

set_option linter.unusedSimpArgs false

namespace NLogN

open IntMul

variable {P : ℕ} {I : List Sym} {iprog : Fin P → IStmt HS P}

namespace RS

/-! ## Frame -/

/-! ## Compilation -/

/-! ## Sequences and reasoning lemmas -/

/-! ## Cost bounds -/

theorem asgCost_le (i : Fin HS) (e : Ex HS) (s : IState HS) (B : ℕ) (h : e.Bnd s B)
    (hi : zsz (s.top i) ≤ B) : asgCost i e s ≤ (e.size + 1) * 10 ^ 7 * (B + 2) ^ 2 := by
  have h1 := Ex.cost_le_r6c s B e h
  have h2 := Ex.zsz_eval_le s B e h
  simp only [asgCost, moveCost]; nlinarith

end RS

end NLogN

end

section
-- Source module: NLogN.HExprDer
/-!
# Derived expressions

Floor division and modulus for an arbitrary numerator and a positive divisor, built from `sel`
and `ndiv`/`nmod`.
-/

namespace NLogN

open IntMul

namespace Ex

theorem floor_div_neg (a b : ℤ) (hb : 0 < b) : a / b = -((-a + (b + -1)) / b) := by
  have h1 := Int.ediv_mul_le (-a + (b + -1)) hb.ne'
  have h2 := Int.lt_ediv_add_one_mul_self (-a + (b + -1)) hb
  set q := (-a + (b + -1)) / b
  have := (Int.ediv_emod_unique'' (a := a) (b := b) (r := a + b * q) (q := -q) hb.ne').2
    ⟨by ring, by nlinarith, by rw [abs_of_pos hb]; nlinarith⟩
  exact this.1

theorem fdiv_eval (s : IState HS) (a b : Ex HS) (hb : 0 < b.eval s) :
    (fdiv a b).eval s = a.eval s / b.eval s := by
  simp only [fdiv, eval]
  split_ifs with h
  · exact (floor_div_neg _ _ hb).symm
  · rfl

theorem fdiv_ok (s : IState HS) (a b : Ex HS) (ha : a.ok s) (hb' : b.ok s) (hb : 0 < b.eval s) :
    (fdiv a b).ok s := by
  simp only [fdiv, ok, eval]
  refine ⟨ha, ?_⟩
  split_ifs with h
  · exact ⟨⟨ha, hb', trivial⟩, hb', by omega, hb⟩
  · exact ⟨ha, hb', by omega, hb⟩

end Ex

end NLogN

end

section
-- Source module: NLogN.HNumCall
/-!
# Tools for composing numeric register programs

* general integer division `Ex.ediv`, remainder `Ex.emod`, and `Ex.max0` (`toNat`),
* "top-based" wrappers of `piP`, `log2P`, `expP` usable as conditional `simp` lemmas,
* crude size bounds on `piBits`, `log2Bits`, `expBits`.
-/

set_option linter.unusedSimpArgs false

namespace NLogN

open IntMul RS IntMul.HvdH

namespace Ex

@[simp] theorem ediv_eval (s : IState HS) (a b : Ex HS) :
    (ediv a b).eval s = a.eval s / b.eval s := by
  rcases lt_trichotomy (b.eval s) 0 with h | h | h
  · have h1 : ¬ (-b.eval s < 0) := by omega
    simp only [ediv, eval, h1, h, if_false, if_true]
    rw [fdiv_eval s a _ (by simp only [eval]; linarith)]
    simp only [eval, Int.ediv_neg, neg_neg]
  · simp [ediv, eval, h]
  · have h1 : -b.eval s < 0 := by omega
    simp only [ediv, eval, h1, if_true]
    exact fdiv_eval s a b h

@[simp] theorem emod_eval (s : IState HS) (a b : Ex HS) :
    (emod a b).eval s = a.eval s % b.eval s := by
  simp only [emod, eval, ediv_eval]
  rw [Int.emod_def]; ring

theorem ediv_ok {s : IState HS} {a b : Ex HS} (ha : a.ok s) (hb : b.ok s) : (ediv a b).ok s := by
  rcases lt_trichotomy (b.eval s) 0 with h | h | h
  · have h1 : ¬ (-b.eval s < 0) := by omega
    simp only [ediv, ok, eval, h1, h, if_false, if_true]
    exact ⟨hb, hb, fdiv_ok s a _ ha hb (by simp only [eval]; linarith)⟩
  · simp [ediv, ok, eval, h, hb]
  · have h1 : -b.eval s < 0 := by omega
    simp only [ediv, ok, eval, h1, if_true]
    exact ⟨hb, fdiv_ok s a b ha hb h⟩

end Ex

/-! ## Top-based wrappers -/

/-! ## Size bounds of the constants -/

end NLogN

end

section
-- Source module: NLogN.SSDif
/-!
# The forward FFT (decimation in frequency) over `ℤ/(2^n+1)`, as a recursive procedure
-/

namespace NLogN

open IntMul

/-! ## Reference functions -/

/-! ## The butterfly step -/

variable {P : ℕ} {I : List Sym} {iprog : Fin P → IStmt NS P}

end NLogN

namespace NLogN

open IntMul

variable {P : ℕ} {I : List Sym} {iprog : Fin P → IStmt NS P}

theorem zsz_natCast_le_self (k : ℕ) : zsz (k : ℤ) ≤ k := by
  apply zsz_le_of_lt; simp only [Int.natAbs_natCast]; exact Nat.lt_two_pow_self

end NLogN

end

section
-- Source module: Solutions.ButterflyGuard
/-!
# §5.4 An explicit guard-width bound (and the numeric parts of Proposition 14)

1. The constants `E = 64(W+m+1)³`, `B = s+E`, `ν = min{j ≥ 1 : m^j ≥ B}`, `C1 = max{20, ν+3}`,
   `C0 = 128 m B²` for the complex network `m = 10⁶`, `W_c`, `s_c`: `m^10 < B < m^11`, `ν = 11`,
   `C1 = 20`.
2. The depth recurrence `A(e) ≤ s A(e/m) + E`, `A ≤ 8e` at leaves, gives `A ≤ 8 d B^{k+1}
   ≤ 8 B² d^{ν+1}`, and the final chain `8mB²d^{ν+2}+8d+t+9 ≤ … < 128mB²d^{ν+3} ≤ C0 d^{C1}`.
3. Dyadic growth: a value at arithmetic depth `ℓ` (add, sub, halve, negate, multiply by `i`)
   lies in `2^{-(p+ℓ)} ℤ[i]` with modulus `≤ 2^ℓ`; hence `p+Δ` fractional and `Δ+3` integer and
   sign bits represent it exactly.
4. The final truncation `Q_p` toward zero: `|Q_p z| ≤ |z|`, `|Q_p z - z| < √2 · 2^{-p}`.
5. The width bookkeeping of the proof of Proposition 14.
-/

namespace Butterfly

open Complex

/-! ## 1. The constants -/

/-- `E = 64 (W + m + 1)³`. -/
def guardE (m W : ℕ) : ℕ := 64 * (W + m + 1) ^ 3
/-- `B = s + E`. -/
def guardB (m W s : ℕ) : ℕ := s + guardE m W

/-! ## 2. Depth recurrence and the chain -/

/-! ## 3. Dyadic depth invariant -/

/-! ## 4. Final truncation toward zero -/

/-! ## 5. Width bookkeeping in the proof of Proposition 14 -/

end Butterfly

end

section
-- Source module: Solutions.KFMBaseSap
open Complex

namespace IntMul.HvdH.SPf
open IntMul.HvdH

/-- The exact weight `e_j = e^{-π (j - sk/t)² / α²}`. -/
noncomputable def ew (s t α k : ℕ) (j : ℤ) : ℝ := Real.exp (-(Real.pi * (j - (s * k : ℝ) / t) ^ 2 / α ^ 2))

end IntMul.HvdH.SPf

end

section
-- Source module: Solutions.J6AMath1
/-!
# A lines (PR #5 block chirp), part 1: the block-chirp `S'` map and its error

The output axis (length `t`) is cut into blocks of `B = 2H` consecutive indices. For an output `k`
the block centre is `k_c = ⌊k/B⌋ B + H`, the offset `κ = k - k_c ∈ [-H, H)`, and the input centre
is `j_c = round(s k_c / t)`. With `A = t j - s k_c` the exponent of the Gaussian weight splits as

`(j - s k/t)² = ((t-s) A² + s (A - t κ)²) / t³ - s (t-s) κ² / t²`   (`expo_split`)

so `e^{-π(j - sk/t)²/α²} = a · g · C` with an input scaling `a` (per block), a correlation kernel
`g` (per block, a function of `j - κ`), and an output scaling `C = e^{π s(t-s)κ²/(t²α²)}` that is the
same for every block (`ew_split`). The machine computes, for every output,

`chirpA U k = (Σ_{|j - j_c| < W} aT_j gT_{j-κ} U_{j mod s}) · cT_κ  tdiv  2^{3P-σ}`

with `aT ≈ 2^P a`, `gT ≈ 2^P g` (`expNeg`) and `cT ≈ 2^{P-σ} C / (2α)` (`expPos`, `betaT`).

**`chirpA_err`**: for every input in the disc and every `k`,
`‖chirpA U k - 2^p S'(U/2^p)_k‖ ≤ 12 W 2^{p+σ} / 2^P + 2`, as soon as the window contains the
`Sapprox` window (`H + m + 3 ≤ W`) and `e^{π s(t-s)H²/(t²α²)} ≤ 2^σ`.
-/

set_option linter.unusedVariables false
set_option linter.unusedSectionVars false

open Complex

namespace J6A

open IntMul.HvdH IntMul.HvdH.SPf

/-! ## 1. Block geometry and tables -/

/-- The centre of the block of `k` (blocks of length `2H`). -/
def kcB (H k : ℕ) : ℕ := k / (2 * H) * (2 * H) + H

/-- The offset `κ = k - k_c` of `k` in its block. -/
def kapB (H k : ℕ) : ℤ := (k : ℤ) - kcB H k

/-- The input centre `round(s k_c / t)` of the block of `k`. -/
def jcB (s t H k : ℕ) : ℤ := (((2 * s * kcB H k + t) / (2 * t) : ℕ) : ℤ)

/-- Numerator of the input-scaling exponent: `(t-s) (t j - s k_c)²` (denominator `t³α²`). -/
def aNum (s t H k : ℕ) (j : ℤ) : ℕ := (t - s) * ((t : ℤ) * j - s * kcB H k).natAbs ^ 2

/-- Numerator of the kernel exponent: `s (t (j - κ) - s k_c)²` (denominator `t³α²`). -/
def gNum (s t H k : ℕ) (j : ℤ) : ℕ := s * ((t : ℤ) * (j - kapB H k) - s * kcB H k).natAbs ^ 2

/-- Numerator of the output-scaling exponent: `s (t-s) κ²` (denominator `t²α²`). -/
def cNum (s t H k : ℕ) : ℕ := s * (t - s) * (kapB H k).natAbs ^ 2

/-- `≈ 2^P e^{-π (t-s)(tj - s k_c)²/(t³α²)}`. -/
def aT (P s t α H k : ℕ) (j : ℤ) : ℤ := expNeg P (t ^ 3 * α ^ 2) (aNum s t H k j)

/-- `≈ 2^P e^{-π s(t(j-κ) - s k_c)²/(t³α²)}`. -/
def gT (P s t α H k : ℕ) (j : ℤ) : ℤ := expNeg P (t ^ 3 * α ^ 2) (gNum s t H k j)

/-- `≈ 2^{P-σ} e^{π s(t-s)κ²/(t²α²)} / (2α)`. -/
def cT (P σ s t α H k : ℕ) : ℤ :=
  (expPos P (t ^ 2 * α ^ 2) (cNum s t H k) σ * betaT P α).tdiv (2 ^ P)

/-- One term `aT_j gT_{j-κ} U_{j mod s}` (before the output scaling). -/
def cTerm (P s t α H k : ℕ) (U : ℕ → ℤ × ℤ) (j : ℤ) : ℤ × ℤ :=
  let w := aT P s t α H k j * gT P s t α H k j
  let V := U (j % s).toNat
  (w * V.1, w * V.2)

/-- **The block-chirp `S'` map** (window `|j - j_c| < W`, precision `P`, scaling exponent `σ`). -/
def chirpA (P σ W s t α H : ℕ) (U : ℕ → ℤ × ℤ) (k : ℕ) : ℤ × ℤ :=
  let c := cT P σ s t α H k
  let jc := jcB s t H k
  (J6Seg.fmul (3 * P - σ) c (∑ i ∈ Finset.range (2 * W), (cTerm P s t α H k U (jc - W + i)).1),
    J6Seg.fmul (3 * P - σ) c (∑ i ∈ Finset.range (2 * W), (cTerm P s t α H k U (jc - W + i)).2))

/-! ## 2. Geometry facts -/

/-! ## 3. The exponent split -/

/-! ## 4. One term -/

end J6A

end

section
-- Source module: Solutions.KFMBaseEap
open Complex

namespace IntMul.HvdH.EPf
open IntMul.HvdH

/-- The exact weight `e_j = e^{-π (j - sk/t)² / α²}`. -/
lemma hasSum_geom_int {r : ℝ} (hr0 : 0 ≤ r) (hr1 : r < 1) :
    HasSum (fun n : ℤ => r ^ n.natAbs) ((1 + r) / (1 - r)) := by
  have h1 : HasSum (fun n : ℕ => r ^ ((n : ℤ).natAbs)) (1 / (1 - r)) := by
    simpa [one_div] using hasSum_geometric_of_lt_one hr0 hr1
  have h2 : HasSum (fun n : ℕ => r ^ ((-((n : ℤ) + 1)).natAbs)) (r / (1 - r)) := by
    have e : (fun n : ℕ => r ^ ((-((n : ℤ) + 1)).natAbs)) = fun n : ℕ => r * r ^ n := by
      funext n; rw [show (-((n : ℤ) + 1)).natAbs = n + 1 by omega, pow_succ, mul_comm]
    rw [e, div_eq_mul_inv]; exact (hasSum_geometric_of_lt_one hr0 hr1).mul_left r
  have : 1 - r ≠ 0 := by linarith
  rw [show (1 + r) / (1 - r) = 1 / (1 - r) + r / (1 - r) by field_simp]
  exact h1.of_nat_of_neg_add_one h2

lemma exp_neg_pi_half_le : Real.exp (-Real.pi / 2) ≤ 1 / 4.3 := by
  have hpi := Real.pi_gt_d2
  have h := Real.quadratic_le_exp_of_nonneg (x := Real.pi / 4) (by positivity)
  have e : Real.exp (Real.pi / 2) = Real.exp (Real.pi / 4) * Real.exp (Real.pi / 4) := by
    rw [← Real.exp_add]; ring_nf
  have hq : 2.09 ≤ Real.exp (Real.pi / 4) := by nlinarith
  have h2 : 4.3 ≤ Real.exp (Real.pi / 2) := by rw [e]; nlinarith
  rw [show -Real.pi / 2 = -(Real.pi / 2) by ring, Real.exp_neg, inv_eq_one_div,
    div_le_div_iff₀ (Real.exp_pos _) (by norm_num)]
  linarith

lemma bOf_bounds (s t : ℕ) (hs : 0 < s) (l : ℤ) :
    -(s : ℤ) ≤ 2 * bOf s t l ∧ 2 * bOf s t l < s := by
  have hs2 : (0 : ℤ) < 2 * s := by omega
  have h1 := Int.ediv_mul_le (2 * t * l + s) hs2.ne'
  have h2 := Int.lt_ediv_add_one_mul_self (2 * t * l + s) hs2
  unfold bOf
  generalize (2 * (t : ℤ) * l + s) / (2 * s) = q at h1 h2 ⊢
  constructor <;> nlinarith

lemma eEx_nonneg (s t α : ℕ) (l h : ℤ) : 0 ≤ eEx s t α l h := by
  unfold eEx; split_ifs
  · exact le_rfl
  · exact (Real.exp_pos _).le

end IntMul.HvdH.EPf

end

section
-- Source module: Solutions.KFMBaseJap
open Complex Real

namespace IntMul.HvdH.JPf
open IntMul.HvdH

lemma hasSum_geom_int {r : ℝ} (hr0 : 0 ≤ r) (hr1 : r < 1) :
    HasSum (fun n : ℤ => r ^ n.natAbs) ((1 + r) / (1 - r)) := by
  have h1 : HasSum (fun n : ℕ => r ^ ((n : ℤ).natAbs)) (1 / (1 - r)) := by
    simpa [one_div] using hasSum_geometric_of_lt_one hr0 hr1
  have h2 : HasSum (fun n : ℕ => r ^ ((-((n : ℤ) + 1)).natAbs)) (r / (1 - r)) := by
    have e : (fun n : ℕ => r ^ ((-((n : ℤ) + 1)).natAbs)) = fun n : ℕ => r * r ^ n := by
      funext n; rw [show (-((n : ℤ) + 1)).natAbs = n + 1 by omega, pow_succ, mul_comm]
    rw [e, div_eq_mul_inv]; exact (hasSum_geometric_of_lt_one hr0 hr1).mul_left r
  have : 1 - r ≠ 0 := by linarith
  rw [show (1 + r) / (1 - r) = 1 / (1 - r) + r / (1 - r) by field_simp]
  exact h1.of_nat_of_neg_add_one h2

/-- Variant excluding `n = 0`: if `f 0 = 0` and `0 ≤ f n ≤ C r^{|n|}` then `∑ f ≤ 2 C r/(1-r)`. -/
lemma tsum_int_le_geom' (f : ℤ → ℝ) (hf0 : ∀ n, 0 ≤ f n) (hz : f 0 = 0) {C r : ℝ}
    (hr0 : 0 ≤ r) (hr1 : r < 1) (hf : ∀ n, n ≠ 0 → f n ≤ C * r ^ n.natAbs) :
    Summable f ∧ ∑' n, f n ≤ 2 * C * r / (1 - r) := by
  have hC : ∀ n : ℤ, n ≠ 0 → 0 ≤ C * r ^ n.natAbs := fun n hn => (hf0 n).trans (hf n hn)
  set g : ℤ → ℝ := fun n => if n = 0 then 0 else C * r ^ n.natAbs
  have hg' := (hasSum_geom_int hr0 hr1).mul_left C
  have hg : HasSum g (C * ((1 + r) / (1 - r)) - C) := by
    have := hg'.update 0 0
    have hgu : g = Function.update (fun i : ℤ => C * r ^ i.natAbs) 0 0 := by
      funext n; by_cases h : n = 0 <;> simp [g, h]
    have e : C * ((1 + r) / (1 - r)) - C = 0 - C * r ^ Int.natAbs 0 + C * ((1 + r) / (1 - r)) := by
      simp only [Int.natAbs_zero, pow_zero, mul_one]; ring
    rw [hgu, e]; exact this
  have hfg : ∀ n, f n ≤ g n := by
    intro n; simp only [g]; split_ifs with h
    · rw [h, hz]
    · exact hf n h
  have hs : Summable f := Summable.of_nonneg_of_le hf0 hfg hg.summable
  refine ⟨hs, (hs.tsum_le_tsum hfg hg.summable).trans (le_of_eq ?_)⟩
  rw [hg.tsum_eq]
  have : 1 - r ≠ 0 := by linarith
  field_simp; ring

/-- Key inequality of Lemma 4.6. -/
lemma key_ineq {q b1 b2 : ℝ} (hq : 1 ≤ q) (hb1 : |b1| ≤ 1 / 2) (hb2 : |b2| ≤ 1 / 2) {h : ℤ}
    (hh : h ≠ 0) : 2 * (q - 1) * (|(h : ℝ)| - 1 / 2) ^ 2 ≤ (q * h + b1) ^ 2 - b2 ^ 2 := by
  have ha : 1 ≤ |(h : ℝ)| := by
    rw [← Int.cast_abs]; exact_mod_cast Int.one_le_abs hh
  have hb2' : b2 ^ 2 ≤ 1 / 4 := by
    have := sq_abs b2 ▸ pow_le_pow_left₀ (abs_nonneg b2) hb2 2; nlinarith
  -- |q h + b1| ≥ q (|h| - 1/2)
  have h1 : q * (|(h : ℝ)| - 1 / 2) ≤ |q * h + b1| := by
    have e1 : |q * (h : ℝ)| = q * |(h : ℝ)| := by rw [abs_mul, abs_of_pos (by linarith)]
    have := abs_sub_abs_le_abs_sub (q * (h : ℝ)) (-b1)
    rw [sub_neg_eq_add, abs_neg, e1] at this
    nlinarith
  have h2 : 0 ≤ q * (|(h : ℝ)| - 1 / 2) := by nlinarith
  have h3 : (q * (|(h : ℝ)| - 1 / 2)) ^ 2 ≤ (q * h + b1) ^ 2 := by
    rw [← sq_abs (q * (h : ℝ) + b1)]; exact pow_le_pow_left₀ h2 h1 2
  rw [mul_pow] at h3
  set X := (|(h : ℝ)| - 1 / 2) ^ 2
  have e1 : 2 * q - 1 ≤ q ^ 2 := by nlinarith [sq_nonneg (q - 1)]
  have e2 : 1 / 4 ≤ X := by simp only [X]; nlinarith
  have e3 : (2 * q - 1) * X ≤ q ^ 2 * X := mul_le_mul_of_nonneg_right e1 (sq_nonneg _)
  nlinarith

lemma exp_neg_two_pi_le : Real.exp (-2 * π) ≤ 1 / 26 := by
  have hpi := Real.pi_gt_d2
  have h := Real.quadratic_le_exp_of_nonneg (x := 2 * π) (by positivity)
  rw [show -2 * π = -(2 * π) by ring, Real.exp_neg]
  rw [inv_eq_one_div, div_le_div_iff₀ (Real.exp_pos _) (by norm_num)]
  nlinarith

lemma exp_neg_pi_half_le : Real.exp (-π / 2) ≤ 1 / 4.3 := by
  have hpi := Real.pi_gt_d2
  have h := Real.quadratic_le_exp_of_nonneg (x := π / 4) (by positivity)
  have e : Real.exp (π / 2) = Real.exp (π / 4) * Real.exp (π / 4) := by
    rw [← Real.exp_add]; ring_nf
  have hq : 2.09 ≤ Real.exp (π / 4) := by nlinarith
  have h2 : 4.3 ≤ Real.exp (π / 2) := by rw [e]; nlinarith
  rw [show -π / 2 = -(π / 2) by ring, Real.exp_neg, inv_eq_one_div,
    div_le_div_iff₀ (Real.exp_pos _) (by norm_num)]
  linarith

/-- Numerical inequality behind Lemma 4.6: for `x ≥ 1`,
`2 e^{-πx/2} / (1 - e^{-2πx}) < 2^{-x}`. -/
lemma numeric_E {x : ℝ} (hx : 1 ≤ x) :
    2 * Real.exp (-π * x / 2) / (1 - Real.exp (-2 * π * x)) < (2 : ℝ) ^ (-x) := by
  have hpi := Real.pi_gt_d2
  have hl := Real.log_two_lt_d9
  have hr : Real.exp (-2 * π * x) ≤ 1 / 26 := by
    refine le_trans (Real.exp_le_exp.mpr ?_) exp_neg_two_pi_le
    nlinarith [Real.pi_pos]
  have hr0 := Real.exp_pos (-2 * π * x)
  set A := Real.exp (-x * Real.log 2) with hA
  have hApos : 0 < A := Real.exp_pos _
  have h2x : (2 : ℝ) ^ (-x) = A := by
    rw [Real.rpow_def_of_pos two_pos, hA]; ring_nf
  rw [h2x, div_lt_iff₀ (by linarith)]
  have hsplit : Real.exp (-π * x / 2) = A * Real.exp (-x * (π / 2 - Real.log 2)) := by
    rw [hA, ← Real.exp_add]; ring_nf
  have hmono : Real.exp (-x * (π / 2 - Real.log 2)) ≤ Real.exp (-(π / 2 - Real.log 2)) := by
    apply Real.exp_le_exp.mpr; nlinarith
  have hval : Real.exp (-(π / 2 - Real.log 2)) = 2 * Real.exp (-π / 2) := by
    rw [show -(π / 2 - Real.log 2) = Real.log 2 + -π / 2 by ring, Real.exp_add,
      Real.exp_log two_pos]
  have hb := exp_neg_pi_half_le
  have : Real.exp (-x * (π / 2 - Real.log 2)) ≤ 2 / 4.3 := by
    rw [hval] at hmono; linarith
  rw [hsplit]
  nlinarith

variable (s t α : ℕ)

lemma bOf_bounds (hs : 0 < s) (l : ℤ) :
    -(s : ℤ) ≤ 2 * bOf s t l ∧ 2 * bOf s t l < s := by
  have hs2 : (0 : ℤ) < 2 * s := by omega
  have h1 := Int.ediv_mul_le (2 * t * l + s) hs2.ne'
  have h2 := Int.lt_ediv_add_one_mul_self (2 * t * l + s) hs2
  unfold bOf
  generalize (2 * (t : ℤ) * l + s) / (2 * s) = q at h1 h2 ⊢
  constructor <;> nlinarith

lemma abs_b_le (hs : 0 < s) (l : ℤ) : |(bOf s t l : ℝ) / s| ≤ 1 / 2 := by
  obtain ⟨h1, h2⟩ := bOf_bounds s t hs l
  have hsR : (0 : ℝ) < s := by exact_mod_cast hs
  have h1' : -(s : ℝ) ≤ 2 * bOf s t l := by exact_mod_cast h1
  have h2' : 2 * (bOf s t l : ℝ) ≤ s := by exact_mod_cast h2.le
  rw [abs_div, abs_of_pos hsR, div_le_iff₀ hsR, abs_le]
  constructor <;> linarith

lemma eEx_nonneg (l h : ℤ) : 0 ≤ eEx s t α l h := by
  unfold eEx; split_ifs
  · exact le_rfl
  · exact (Real.exp_pos _).le

lemma eEx_le_gauss (hs : 0 < s) (hst : s ≤ t) (l h : ℤ) (hh : h ≠ 0) :
    eEx s t α l h ≤ Real.exp (-2 * π * xθ s t α * (|(h : ℝ)| - 1 / 2) ^ 2) := by
  unfold eEx; rw [if_neg hh]
  apply Real.exp_le_exp.2
  have hsR : (0 : ℝ) < s := by exact_mod_cast hs
  have hq : 1 ≤ (t : ℝ) / s := by rw [one_le_div hsR]; exact_mod_cast hst
  have key := key_ineq hq (abs_b_le s t hs l) (abs_b_le s t hs (l + h)) hh
  have hK : (eKz s t l h : ℝ) / s ^ 2 =
      ((t : ℝ) / s * h + (bOf s t l : ℝ) / s) ^ 2 - ((bOf s t (l + h) : ℝ) / s) ^ 2 := by
    unfold eKz; push_cast; field_simp
  have hx : xθ s t α = α ^ 2 * ((t : ℝ) / s - 1) := by
    unfold xθ; rw [Nat.cast_mul, Nat.cast_sub hst]; push_cast; field_simp
  have hpa : (0 : ℝ) ≤ π * α ^ 2 := by positivity
  have := mul_le_mul_of_nonneg_left key hpa
  rw [← hK] at this
  rw [hx]
  have e : π * α ^ 2 * (eKz s t l h : ℝ) / s ^ 2 = π * α ^ 2 * ((eKz s t l h : ℝ) / s ^ 2) := by ring
  rw [e]; nlinarith

lemma eEx_sum (hs : 0 < s) (hst : s ≤ t) (hx : 1 ≤ xθ s t α) (l : ℤ) :
    Summable (eEx s t α l) ∧ ∑' h, eEx s t α l h ≤ (2 : ℝ) ^ (-xθ s t α) := by
  set x := xθ s t α
  have hx0 : 0 ≤ x := by linarith
  have hgeom := tsum_int_le_geom' (eEx s t α l) (eEx_nonneg s t α l) (by simp [eEx])
    (C := Real.exp (3 * π * x / 2)) (r := Real.exp (-2 * π * x)) (Real.exp_pos _).le
    (Real.exp_lt_one_iff.mpr (by nlinarith [Real.pi_pos])) (by
      intro h hh
      refine (eEx_le_gauss s t α hs hst l h hh).trans ?_
      rw [← Real.exp_nat_mul, ← Real.exp_add]
      apply Real.exp_le_exp.mpr
      have ha : (h.natAbs : ℝ) = |(h : ℝ)| := by rw [Nat.cast_natAbs, Int.cast_abs]
      rw [ha]
      nlinarith [mul_nonneg (mul_nonneg Real.pi_pos.le hx0) (sq_nonneg (|(h : ℝ)| - 1))])
  refine ⟨hgeom.1, hgeom.2.trans (le_of_lt ?_)⟩
  have e : 2 * Real.exp (3 * π * x / 2) * Real.exp (-2 * π * x) = 2 * Real.exp (-π * x / 2) := by
    rw [mul_assoc, ← Real.exp_add]; ring_nf
  rw [e]; exact numeric_E hx

end IntMul.HvdH.JPf

end

section
-- Source module: Solutions.ResampleHvdH
/-! Copy (renamed namespace `IntMul.OpenAI.RH`) of the base thread's HvdH Prop 4.7(i)
(`intmul/base/Sol_prop_4_7_i.lean`): Gaussian resampling identity, maps S, T, C, D, P_s, P_t. -/

set_option linter.unusedSectionVars false
set_option linter.unusedSimpArgs false

section Resample

open Complex Real

namespace IntMul.OpenAI.RH

/-! ### Definitions (HvdH §2.4, §4.1) -/

/-! ### Gaussian sums -/

/-! ### Characters -/

/-! ### Linearity -/

/-! ### The resampling identity (HvdH Theorem 4.2) -/

end IntMul.OpenAI.RH

end Resample

section Bounds

open Real

namespace IntMul.OpenAI.RH

end IntMul.OpenAI.RH

end Bounds

section Maps

open Complex Real

namespace IntMul.OpenAI.RH

section maps
variable {s t : ℕ} [NeZero s] [NeZero t]

variable (s t) in
/-- Row selection `C : ℂ^t → ℂ^s`, `(C v)_l = v_{[t l / s]}`. -/
noncomputable def rowSel (v : ZMod t → ℂ) : ZMod s → ℂ :=
  fun l => v ((rnd ((t : ℝ) * l.val / s) : ℤ) : ZMod t)

end maps

end IntMul.OpenAI.RH

end Maps

section Prop47

open Complex Real

namespace IntMul.OpenAI.RH

end IntMul.OpenAI.RH

end Prop47

end

section
-- Source module: NLogN.HSort
/-!
# Stable radix sort of fixed-width records on a register stack

* `mvK X Y C n`: move `n` integers from `X` to `Y` (reversed); `mvK_spec`, `mvK_cost`.
* `radixPass k`, `radixSort kb`: the pure stable partition pass / LSD radix sort on the first
  component (the key) of each record; `radixSort_keys`: for keys forming a permutation of
  `0, …, N-1` (with `N ≤ 2^kb`) the record with key `q` ends up at position `q`.
* `passP w`, `sortP w`: the register programs (records of width `w ≥ 1` on stack `0`, record
  count on `1`, key bits `kb` on `2`; scratch `46 … 53`); `sortP_spec`, `sortP_cost`.
All per-element code uses `Ex.Lin` expressions only; costs are bounded with the linear lemmas.
-/

set_option linter.unusedSimpArgs false

namespace NLogN

open IntMul RS

theorem mvK_cost {X Y C : Fin HS} (hXY : X ≠ Y) (hXC : X ≠ C) (hYC : Y ≠ C) (n : Ex HS)
    (s : IState HS) (L : List ℤ) (R : List Item) (hX : s.G X = ints L ++ R)
    (hn : n.eval s = L.length) (B : ℕ) (hL : Bnd B L) (hnl : n.Lin) (hnb : n.Bnd s B) :
    (mvK X Y C n).cost s ≤ (n.size + 1) * 2000 * (B + 2) + L.length * (14000 * (B + 2)) + 32 := by
  have H := mvK_states hXY hXC hYC n s L R hX hn
  have := H.cost_le_lin B (4000 * (B + 2) + 2 * B + 2) hnl hnb ?_
  · rw [hn, Int.toNat_natCast] at this
    have h2 : L.length * (4000 * (B + 2) + 2 * B + 2 + 8000 * (B + 2) + 30) ≤
        L.length * (14000 * (B + 2)) := Nat.mul_le_mul_left _ (by omega)
    exact le_trans this (by omega)
  intro k hk
  rw [hn, Int.toNat_natCast] at hk
  have hd : L.drop k = L[k] :: L.drop (k + 1) := List.drop_eq_getElem_cons hk
  have e1 : (withCtr C s (n.eval s).toNat k (mvSt_s X Y s L R k)).G X =
      .int L[k] :: (ints (L.drop (k + 1)) ++ R) := by
    simp only [withCtr, mvSt_s, IState.upd_G, if_neg hXY, if_neg hXC, ite_true]
    rw [hd]; rfl
  have hb : zsz L[k] ≤ B := hL _ (List.getElem_mem hk)
  have h1 := pushCost_le_lin (.reg X) (withCtr C s (n.eval s).toNat k (mvSt_s X Y s L R k)) B
    (by simp [Ex.Lin]) (by simpa [Ex.Bnd, IState.top_eq_r6c e1] using hb)
  simp only [RS.cost, run_seq]
  have e2 : ((push Y (.reg X)).run (withCtr C s (n.eval s).toNat k (mvSt_s X Y s L R k))).top X =
      L[k] := IState.top_eq_r6c (L := ints (L.drop (k + 1)) ++ R)
        (by simp only [run_push_G, if_neg hXY, e1])
  rw [e2]
  have : (Ex.reg (S := HS) X).size = 1 := rfl
  rw [this] at h1
  omega

/-! ## Radix sort: pure functions -/

/-! ## One pass as a register program

Registers: data `0`, record count `1`; scratch `46` (ones), `47` (zeros), `48`, `49` (their
element counts), `50` (loop counter), `51` (counter of `mvK`), `52` (holds `2^k`). -/

/-! ## The radix sort program

Data on `0`, record count on `1`, number of key bits on `2`; scratch `46 … 53`. -/

/-! ## Costs -/

/-! ## Correctness of radix sort for distinct keys -/

end NLogN

end

section
-- Source module: Solutions.KILBase
/-!
# KIL base: register programs on "managed" stacks

`encK s f` is the state `s` with the list of integers `f j` pushed on top of every stack `j`
(`f j = []` for unmanaged stacks).  All programs of the KIL files are specified on such states;
updates are `Function.update`s of `f`, which keeps the bookkeeping small.

`Sp p s t c` : the register program `p` is executable from `s` (`RS.ok`), ends in `t`, and costs
at most `c`.  `Sp.irun` turns it into an `IRuns` of the compiled item-language code.
-/

set_option linter.unusedSimpArgs false

namespace NLogN

open IntMul RS

namespace KIL

/-! ## Elementary statements -/

/-! ## Moving integers between stacks -/

theorem Sp.mvK_r6c {X Y C : Fin HS} (hXY : X ≠ Y) (hXC : X ≠ C) (hYC : Y ≠ C) (hX5 : X.val < 56)
    (hY5 : Y.val < 56) (hC5 : C.val < 56) (s : IState HS) (f : Fin HS → List ℤ) (n : Ex HS)
    (hnc : n.Clean) (hno : n.ok (encK s f)) (L l : List ℤ) (hX : f X = L ++ l)
    (hn : n.eval (encK s f) = L.length) (B : ℕ) (hL : Bnd B L) (hnl : n.Lin)
    (hnb : n.Bnd (encK s f) B) :
    Sp (mvK X Y C n) (encK s f)
      (encK s (Function.update (Function.update f X l) Y (L.reverse ++ f Y)))
      ((n.size + 1) * 2000 * (B + 2) + L.length * (14000 * (B + 2)) + 32) := by
  have hXs : (encK s f).G X = ints L ++ (ints l ++ s.G X) := by simp [hX, ints_append]
  have H := mvK_spec hXY hXC hYC hX5 hY5 hC5 n hnc (encK s f) L _ hXs hno hn
  refine ⟨H.1, ?_, mvK_cost hXY hXC hYC n (encK s f) L _ hXs hn B hL hnl hnb⟩
  rw [H.2]
  refine IState.ext' (fun j => ?_) rfl rfl
  by_cases h1 : j = Y
  · subst h1; simp [Ne.symm hXY, ints_append, ints_reverse]
  by_cases h2 : j = X
  · subst h2; simp [hXY]
  · simp [h1, h2]

/-! ## Counted loops through a family of managed states -/

/-! ## Expression side conditions on managed states -/

end KIL

end NLogN

end

section
-- Source module: Solutions.KIL3Tab
/-!
# KIL3: a table of `W` streams on one stack

The `W` streams of the network are kept on stack `0`, one after the other (stream `0` on top):
`tab W F = F 0 ++ F 1 ++ ⋯ ++ F (W − 1)`, every stream of the same length `v`.

Operations (all `O(W · v)` item moves, plus the cost of the applied code):

* `applyC v w code`: run `code` (any stack-`0` program) on stream `w` in place
  (`applyC_runs`): the streams before `w` are parked on stack `15` and moved back.
* `copyC v w`: push a reversed copy of stream `w` on stack `16` (`copyC_sp`).
* `buildC v l`: push on stack `20` the (reversed) pointwise sum mod 2 of the streams `l`.
* `rebuildC W v T`: replace the table by the streams `w ↦ Σ_{w' ∈ T w} F w'` (`rebuildC_sp`).

Scratch stacks: `15` (park), `16` (copy), `17` (temp), `19` (zip out), `20` (new table),
`21` (loop counter).
-/

set_option linter.unusedSimpArgs false
set_option linter.unnecessarySeqFocus false

namespace NLogN.KIL

open IntMul RS AModel

/-! ## Tables -/

/-! ## Moves with literal counts -/

/-! ## Applying a stream code in place -/

/-! ## Cost unit -/

/-! ## Copying a stream -/

theorem TabB.le {W v Bb : ℕ} (h : TabB W v Bb) {n : ℕ} (hn : n ≤ W * v) : zsz (n : ℤ) ≤ Bb :=
  le_trans (zsz_natCast_le' hn) h.vol

/-! ## Accumulating a stream on stack `20` -/

/-! ## Building a summed stream, rebuilding the table -/

end NLogN.KIL

end

section
-- Source module: Solutions.KIL5Swap
/-!
# KIL5: the level-`0` interchange with run-time sizes, and the descriptor stacks

`prog C q 0 P L A B Cs` is the interchange of the two digits (fields of size `q`) of
`P :: L :: swShape A B Cs q`.  KIL3 realised it by `swapB`, two Lemma 3 field moves whose size
parameters are *literals*.  Here the size parameters are *expressions*, evaluated at run time:

* `withPush l p`: push the expressions `l`, run `p`, pop them (`Sp.withPush`);
* `swapE a q b c`: the two field moves with parameters read from expressions (`swapE_sp`):
  on `A0 ++ q :: (B ++ q :: Cs)` it computes the digit interchange when `a, q, b, c` evaluate to
  `vol A0, q, vol B, vol Cs`;
* the descriptor stacks `dJ … dT` (level, `P`, `L`, `vol A`, `vol B`, `vol Cs`, tower of moduli),
  `DescIn`, and `lev0U`, the uniform level-`0` code reading the descriptor (`lev0U_sp`).
-/

set_option linter.unusedSimpArgs false
set_option linter.unnecessarySeqFocus false

namespace NLogN.KIL

open IntMul RS AModel AModelL12

/-! ## Push, run, pop -/

/-! ## Good expressions -/

theorem RdOK.ne0 {r : Fin HS} (h : RdOK r) : r ≠ 0 := by
  rintro rfl; rcases h with h | h
  · exact absurd h (by decide)
  · exact absurd h.1 (by decide)

/-! ## The two field moves with run-time parameters -/

/-! ## The descriptor stacks -/

/-! ## The uniform level-`0` code -/

end NLogN.KIL

end

section
-- Source module: Solutions.KIL6Pad
/-!
# KIL6: one-field padding and extraction as IL loops

* `padP`: on stack `0` an array `xs` of `X · (F · Y)` entries (`[X, F, Y]`, row-major); it
  produces `padL X F F' Y 0 xs` (every block of `F · Y` entries followed by `(F' - F) · Y`
  zeros).  Per block: move `F · Y` entries to the buffer stack `15` (`mvK`), push the zeros;
  finally move everything back (the two reversals cancel).
* `unpadP`: the inverse scan `unpadL X F F' Y` (move `F · Y`, drop `(F' - F) · Y`).

The four counts are read from the parameter registers `16 … 19` (pushed by the caller); the
loop counters are `20`, `21`.  Both cost `O((V' + 1)(Bb + 2))`, `V'` the padded volume
(`padP_sp`, `unpadP_sp`).
-/

set_option linter.unusedSimpArgs false
set_option linter.unnecessarySeqFocus false

namespace NLogN.KIL

open IntMul RS AModel AModelL12

/-! ## Loops with register counts -/

theorem Sp.zeroR {Y C r : Fin HS} (hYC : Y ≠ C) (hY5 : Y.val < 56) (hC5 : C.val < 56)
    (hr5 : r.val < 56) (s : IState HS) (f : Fin HS → List ℤ) (n : ℕ) (hfC : f C = [])
    {Lr : List Item} (hr : (encK s f).G r = .int (n : ℤ) :: Lr) (Bb : ℕ)
    (hn : zsz (n : ℤ) ≤ Bb) :
    Sp (zeroKE Y C (.reg r)) (encK s f) (encK s (Function.update f Y (List.replicate n 0 ++ f Y)))
      (cZeroR n Bb) := by
  have H := Sp.rep (C := C) (n := .reg r) (s := s)
    (F := fun k => Function.update f Y (List.replicate k 0 ++ f Y)) (N := n) (B := Bb)
    (c := (1 + 1) * 2000 * (Bb + 2)) (b := .push Y (.lit 0)) hC5
    (fun k => by simp [Function.update_apply, hYC.symm, hfC])
    (by simp [RS.writes, hYC.symm]) (clean_reg_r6c hr5)
    (by simp only [List.replicate_zero, List.nil_append, Function.update_eq_self]
        exact IState.isInt_of_r6c hr) trivial
    (by simp only [List.replicate_zero, List.nil_append, Function.update_eq_self, Ex.Bnd]
        rw [IState.top_eq_r6c hr]; exact hn)
    (by simp only [List.replicate_zero, List.nil_append, Function.update_eq_self, Ex.eval]
        rw [IState.top_eq_r6c hr]) ?_
  · have e0 : Function.update f Y (List.replicate 0 0 ++ f Y) = f := by simp
    simp only [e0] at H
    exact H.mono (by simp [cZeroR, Ex.size])
  · intro k hk
    have h1 := Sp.push s (Function.update (Function.update f Y (List.replicate k 0 ++ f Y)) C
      [(n : ℤ) - k]) (i := Y) (e := .lit 0) (B := Bb) hY5 (clean_lit_r6c _) trivial trivial
      (by simp [Ex.Bnd, zsz_zero_r6c])
    refine Sp.mono (Sp.congr h1 ?_) (by simp [Ex.size])
    congr 1
    funext j
    simp only [Function.update_apply, Ex.eval]
    by_cases h1 : j = Y
    · simp only [h1, if_true, if_neg hYC, List.replicate_succ, List.cons_append]
    · simp only [h1, if_false]

/-! ## The padding program -/

end NLogN.KIL

end

section
-- Source module: Solutions.KIL6NE
/-!
# KIL6: natural-number expressions

`NE` is a small language of natural-number expressions over registers (`+`, `-` (truncation
excluded by hypothesis), `*`, `/`, `%`, `if a < b then x else y`), compiled to `Ex HS`.
`NE.good`: on a state whose registers hold the naturals `v i`, if every intermediate value is
`≤ M` and `zsz M ≤ B`, then the compiled expression is clean, well defined, `B`-bounded and
evaluates to the natural value.  This turns every descriptor computation into arithmetic.
-/

set_option linter.unusedSimpArgs false
set_option linter.unusedVariables false

namespace NLogN.KIL

open IntMul RS AModel AModelL12

namespace NE

theorem lin : ∀ {e : NE}, e.IsLin → e.toEx.Lin
  | reg _, _ | lit _, _ => trivial
  | add a b, h => ⟨lin h.1, lin h.2⟩
  | sub a b, h => ⟨lin h.1, lin h.2⟩
  | lt a b x y, h => ⟨⟨lin h.1, lin h.2.1⟩, lin h.2.2.1, lin h.2.2.2⟩

end NE

/-! ## Managed-state statements -/

end NLogN.KIL

end

section
-- Source module: Solutions.KIL6Reg
/-!
# KIL6: register files

`regState F v D`: the managed state `F` with the natural `v i` pushed on every register
`i ∈ D`.  Assignment lists `asgL L` (each an `NE` expression) act on the valuation `v` only
(`Sp.asgs`); counted loops of assignment lists iterate the valuation map (`Sp.repA`).  `pushZs` /
`popDs` create and remove the register file.
-/

set_option linter.unusedSimpArgs false
set_option linter.unusedVariables false

namespace NLogN.KIL

open IntMul RS AModel AModelL12

/-! ## Assignment lists -/

/-! ## Creating and removing the file -/

theorem regState_nil (F : Fin HS → List ℤ) (v : Fin HS → ℕ) : regState F v [] = F := by
  funext j; simp [regState]

theorem Sp.pushZ (s : IState HS) (B : ℕ) :
    ∀ (D : List (Fin HS)) (F : Fin HS → List ℤ), (∀ i ∈ D, i.val < 56) → D.Nodup →
      Sp (pushZs D) (encK s F) (encK s (regState F (fun _ => 0) D)) (D.length * 4000 * (B + 2) + 1)
  | [], F, _, _ => by rw [regState_nil]; exact (Sp.skip _).mono (by simp)
  | i :: D, F, h56, hnd => by
      have h1 := Sp.push s F (i := i) (e := .lit 0) (B := B) (h56 i (by simp)) (clean_lit_r6c 0) trivial
        trivial (by simp [Ex.Bnd, zsz_zero_r6c])
      have h2 := Sp.pushZ s B D (Function.update F i ((Ex.lit (0 : ℤ)).eval (encK s F) :: F i))
        (fun j hj => h56 j (by simp [hj])) (List.nodup_cons.1 hnd).2
      refine ((h1.seq_r6c h2).congr ?_).mono ?_
      · congr 1; funext j
        simp only [regState, Function.update_apply, Ex.eval, List.mem_cons]
        by_cases hj : j = i
        · subst hj; simp [(List.nodup_cons.1 hnd).1]
        · simp [hj]
      · simp only [Ex.size, List.length_cons]; ring_nf; omega

end NLogN.KIL

end

section
-- Source module: Solutions.KFCBase
/-!
# KFC base: fine-cost item-language loops

Item-language (`IStmt`) building blocks whose per-element cost is `O(1 + zsz a)` for the moved
entry `a`, with **no** binary loop counter per element:

* `mvW X Y`   : move integers from `X` to `Y` until the top of `X` is not an integer
                (delimiter-driven, cost `(76 w + 131)` per entry of bit size `≤ w`);
* `uloop C D b`: run `b` once per mark on the *unary counter* `C` (the marks are moved to `D`
                and then moved back, so the counter is restored); overhead `64` per iteration;
* `mvB X Y C D`: move exactly `n` integers, `n` = number of marks on `C`;
* `toU A X K K2`: push `n` marks on `X`, where `n ≥ 0` is the integer on top of `A`
                (binary counter, cost `O((n+1)(zsz n + 2))`; superseded by the linear
                `KFCUnary.toUL`, which the programs use);
* `clrM X`    : drop all marks on top of `X`.
-/

set_option linter.unusedSimpArgs false

namespace NLogN

open IntMul

namespace KFC

variable {Pp : ℕ} {I : List Sym} {iprog : Fin Pp → IStmt HS Pp}

/-! ## Classes and item lists -/

/-! ## `mvW`: move integers until a delimiter -/

/-! ## Unary counted loops -/

/-! ## `mvB`: move a counted block -/

/-! ## `clrM`: drop the marks on top of a stack -/

/-! ## `toU`: binary to unary (prologue only) -/

def isNZ (c : Cls) : Bool := c != .zero

def toU (A X K K2 : Fin HS) : IStmt HS Pp :=
  .seq (.op (.dup A K))
    (.seq (.loop K isNZ (.seq (.op (.mark X)) (.seq (.op (.const K2 (-1))) (.op (.add K K2 K)))))
      (.op (.drop K)))

end KFC

end NLogN

end

section
-- Source module: Solutions.KFEPad
/-!
# KFE: fine padding loops

`padP` / `unpadP` (KIL6Pad) run through `KFD.fcode`: every counted loop is a unary loop, so
the per-element cost is `O(1)` for bit data and the per-block loop overhead is
`O(zsz (F'·Y))`, bounded by the block length.  Result: fine cost
`cPadF V Bb = 2000 (Bb + 2) + 20000 V` (`V` the padded volume) in place of
`cPad V Bb = 10^5 (V+1)(Bb+2)`.

* `FSp.zeroR`, `FSp.dropR`: fine `zeroKE` / `dropKE` with a register count;
* `padP_fsp`, `unpadP_fsp`: the fine specifications (same denotations as `padP_sp`,
  `unpadP_sp`), for bit data (`Bnd 1`);
* `padS_fsp`, `unpadS_fsp`: the stages with their register assignments.
-/

set_option linter.unusedSimpArgs false
set_option linter.unnecessarySeqFocus false
set_option linter.unusedVariables false

namespace NLogN.KFE

open IntMul RS AModel AModelL12 KIL KFD

/-! ## Fine `zeroKE` / `dropKE` with a register count -/

theorem FSp.zeroR {Y C r : Fin HS} (hYC : Y ≠ C) (hY5 : Y.val < 56) (hC5 : C.val < 56)
    (hr5 : r.val < 56) (s : IState HS) (f : Fin HS → List ℤ) (n : ℕ) (hfC : f C = [])
    {Lr : List Item} (hr : (encK s f).G r = .int (n : ℤ) :: Lr) (Bb : ℕ)
    (hn : zsz (n : ℤ) ≤ Bb) :
    FSp (zeroKE Y C (.reg r)) (encK s f) (encK s (Function.update f Y (List.replicate n 0 ++ f Y)))
      (cLp n 8000 Bb) := by
  have e0 : Function.update f Y (List.replicate 0 0 ++ f Y) = f := by simp
  have H := FSp.repReg (C := C) (r := r) (s := s)
    (F := fun k => Function.update f Y (List.replicate k 0 ++ f Y)) (N := n) (Bb := Bb)
    (c := 8000) (b := .push Y (.lit 0)) hC5 hr5
    (fun k => by simp [Function.update_apply, hYC.symm, hfC])
    (by simp [RS.writes, hYC.symm]) (L := Lr) (by simp only [e0]; exact hr) hn ?_
  · simp only [e0] at H
    exact H
  · intro k hk
    have h1 := Sp.push s (Function.update (Function.update f Y (List.replicate k 0 ++ f Y)) C
      [(n : ℤ) - k]) (i := Y) (e := .lit 0) (B := 0) hY5 (clean_lit_r6c _) trivial trivial
      (by simp [Ex.Bnd, zsz_zero_r6c])
    refine FSp.mono (FSp.congr (FSp.of_norep_r6c rfl h1) ?_) (by simp [Ex.size])
    congr 1
    funext j
    simp only [Function.update_apply, Ex.eval]
    by_cases h1 : j = Y
    · simp only [h1, if_true, if_neg hYC, List.replicate_succ, List.cons_append]
    · simp only [h1, if_false]

/-! ## The fine padding loops -/

/-! ## The stages -/

end NLogN.KFE

end

section
-- Source module: Solutions.NetComplexNum
/-!
# The concrete complex network (`h = 100`): neighbour count and the numbers of §3.1, §3.5

* `count_nb_c`: a triple `S ⊂ [n]` has `C(n-3, 3) + 3(n-3)` complex neighbours (distinct triples
  with even intersection): `C(n-3,3)` with intersection `0` and `3(n-3)` with intersection `2`.
* For `T₀` = triples of `Fin 100`: `z_c = 147731`, `|Nb| = 161700 · 147731`, and the wire count
  `W_c = 2N + I(v z_c + c_c) = 1873807244643542670000` (`c_c = 101` central wires `C_i`, `C_*`).
* The numbers of Proposition 8: `L_c = I c_c h`, `s_c = W_c m - 2N + 2L_c =
  1873807244636671267308000000 < W_c m`, `L_c / N = 101/539`, deficit `2N - 2L_c > 0`.
-/

set_option linter.unusedSectionVars false
set_option linter.unusedSimpArgs false
set_option synthInstance.maxSize 512

namespace NetComplex

open Finset

/-! ## 1. Counting complex neighbours -/

/-! ## 2. The concrete network -/

/-- `L_c = I c_c h`. -/
def Lc : ℕ := 3 * 161700 ^ 2 * 101 * 100

end NetComplex

end

section
-- Source module: Solutions.KFLLoop
/-!
# KFL §3: the segment loop (register programs)

All per-record / per-block arithmetic of Lemma 13 is a loop over the `n` segments
`j_i = ρ + iK`: with `u, v` the two chunk values shifted by `2^ρ`, `w = 2^ρ`,

  `segBody K T : acc += w · T(u, v); u /= 2^K; v /= 2^K; w *= 2^K`

repeated `n` times (`segLoop`).  `segLoop_fsp` gives the closed form (for any expression `T` of
the two registers `rU, rV` with semantics `tf`) and the fine cost `O(n (|T| + 1) (Bb+2)^2)`.
Instances: the offset of a line (`offE`), the `T0` correction (`t0E`), the guard test (`grdE`).
-/

set_option linter.unusedSimpArgs false
set_option linter.unusedVariables false

namespace NLogN.KFL

open IntMul RS KIL KFD

/-! ## Registers -/

theorem clean_all (e : Ex HS) (h : ∀ i ∈ e.regs, i.val < 56) : e.Clean := fun i hi hz => by
  have := h i hi
  simp only [ExSt.all, zX, List.mem_cons, List.not_mem_nil, or_false] at hz
  rcases hz with rfl | rfl | rfl | rfl | rfl | rfl | rfl | rfl <;> exact absurd this (by decide)

/-! ## The loop -/
end NLogN.KFL

end

section
-- Source module: Solutions.KFLKey
/-!
# KFL §10: the destination key of a record, in registers

* `segStep`: push `u, v, w, a` on `rU, rV, rW, rA`, run the segment loop, assign a function of
  the accumulator to a target register, pop.  `segStep_fsp` is the generic specification.
* `lineR l`: one packed line (22) on the registers `rX, rY, rZ` (`lineP` on the triple).
* `t0R`: the ideal map `T0` on `rY` (an XOR with `x ∧ selM (f-1)`).
* `keyProg`: from the rank `m` on `rm` push the key `keyR m` on `rK`
  (`x, y, z` by division, the eight lines in `foldr` order, `T0`, and the rank update).
-/

set_option linter.unusedSimpArgs false
set_option linter.unusedVariables false

namespace NLogN.KFL

open IntMul RS KIL KFD Butterfly

/-! ## Expression specifications -/

theorem ESp.congr {e : Ex HS} {t t' : IState HS} {Bb : ℕ} {v : ℤ} (h : ESp e t Bb v)
    (hr : ∀ i ∈ e.regs, t.G i = t'.G i) : ESp e t' Bb v := by
  obtain ⟨h1, h2, h3, _⟩ := Ex.congrG e hr
  exact ⟨h.clean, h2.1 h.ok, (h3 Bb).1 h.bnd, h1 ▸ h.ev⟩

/-! ## The generic segment step -/

/-! ## Expression combinators -/

/-! ## Registers of the key computation -/

/-! ## Size facts -/

/-! ## One packed line on the registers -/

/-! ## The ideal map `T0` on the registers -/

/-! ## The eight lines -/

/-! ## The key of a rank -/

end NLogN.KFL

end

section
-- Source module: Solutions.ResampleParams
/-!
# §9.1–9.2 Fixed rational choice and sizes

* `constants_ok`: the layer and dimension conditions of (38), `margins_ok`: the seven margins
  of (45) with `min gᵢ = 2⁻¹⁸¹ = 2κ`.
* Axis lengths `t_i ∈ {r/2, r}` with `∏ t_i = T`, last axis `r` (`axExp_sum`, `axLen_prod`,
  `axExp_last`), `p/(12d) ≤ ℓ ≤ p/(3d)` (`ell_bounds`), `24n ≤ T p < 48 n` (`Tp_bounds`).
* Source lengths: `(1-2η)t_i < s_i ≤ (1-η)t_i` ⇒ `S > T/2` (`S_gt_half`), `θ_i > η`
  (`theta_gt_eta`), `α⁴θ_i > p` (`alpha4_theta_gt_p`), `γ < 28 d² √b` (`gamma_bound`);
  existence of pairwise coprime odd `s_i` (`sources_exist`, from `RCop.coprime_params`).
* Radix digits (§9.2.3): `2q - 2 < S` (`degree_lt_S`) and `c_j < 2^{3b}` (`coeff_lt`).
-/

open Finset

namespace IntMul.OpenAI.RP

/-! ### (38) and (45) -/

def g7 : ℚ := ε

/-! ### Axis lengths of the padded box -/

lemma ell_spec (L d : ℕ) (hd : 0 < d) : L ≤ d * ell L d ∧ d * ell L d ≤ L + d - 1 := by
  unfold ell
  constructor
  · have := Nat.lt_div_mul_add (a := L + d - 1) hd
    rw [mul_comm]; omega
  · have := Nat.div_mul_le_self (L + d - 1) d
    rw [mul_comm]; exact this

lemma mcnt_lt (L d : ℕ) (hd : 0 < d) : mcnt L d < d := by
  have := ell_spec L d hd; unfold mcnt; omega

lemma axExp_sum (L d : ℕ) (hd : 0 < d) : ∑ i : Fin d, axExp L d i = L := by
  have hs := ell_spec L d hd
  have hm := mcnt_lt L d hd
  have hl : 0 < mcnt L d → 1 ≤ ell L d := by
    intro hpos
    by_contra h; push Not at h
    have : ell L d = 0 := by omega
    unfold mcnt at hpos; rw [this] at hpos; simp at hpos
  have key : ∀ i : Fin d, (axExp L d i : ℤ) =
      (ell L d : ℤ) - (if (i : ℕ) < mcnt L d then 1 else 0) := by
    intro i; unfold axExp
    split_ifs with h
    · have := hl (by omega); push_cast [Nat.cast_sub this]; ring
    · simp
  zify
  rw [Finset.sum_congr rfl (fun i _ => key i), Finset.sum_sub_distrib,
    Fin.sum_univ_eq_sum_range (fun i => if i < mcnt L d then (1 : ℤ) else 0), Finset.sum_boole]
  have hf : (Finset.range d).filter (fun i => i < mcnt L d) = Finset.range (mcnt L d) := by
    ext x; simp only [Finset.mem_filter, Finset.mem_range]; omega
  rw [hf]
  simp only [Finset.sum_const, Finset.card_univ, Fintype.card_fin, Finset.card_range, nsmul_eq_mul]
  unfold mcnt
  push_cast [Nat.cast_sub hs.1]
  ring

/-- `∏ t_i = T` with `t_i = 2^{axExp}`. -/
lemma axLen_prod (L d : ℕ) (hd : 0 < d) : ∏ i : Fin d, 2 ^ axExp L d i = 2 ^ L := by
  rw [Finset.prod_pow_eq_pow_sum, axExp_sum L d hd]

/-- `p/(12d) ≤ ℓ ≤ p/(3d)` for `p = 6b`, `b/2 ≤ L ≤ b`, `1 ≤ d ≤ b`. -/
lemma ell_bounds (L d b : ℕ) (hd : 0 < d) (hdb : d ≤ b) (hL1 : b ≤ 2 * L) (hL2 : L ≤ b) :
    (6 * b : ℝ) / (12 * d) ≤ ell L d ∧ (ell L d : ℝ) ≤ (6 * b : ℝ) / (3 * d) := by
  have hs := ell_spec L d hd
  have hd' : (0 : ℝ) < d := by exact_mod_cast hd
  constructor
  · rw [div_le_iff₀ (by positivity)]
    have : (b : ℝ) ≤ 2 * (d * ell L d) := by
      have : b ≤ 2 * (d * ell L d) := by omega
      exact_mod_cast this
    nlinarith
  · rw [le_div_iff₀ (by positivity)]
    have : ((d * ell L d : ℕ) : ℝ) ≤ 2 * b := by
      have : d * ell L d ≤ 2 * b := by omega
      exact_mod_cast this
    push_cast at this
    nlinarith

/-! ### Source lengths, resampling margins -/

/-- (40): `α ≤ 2(12d²b)^{1/4}` gives `γ = 2dα² < 28 d² √b`. -/
lemma gamma_bound (d b α : ℝ) (hd : 0 < d) (hb : 0 < b) (hα0 : 0 ≤ α)
    (hα : α ≤ 2 * (12 * d ^ 2 * b) ^ ((1 : ℝ) / 4)) :
    2 * d * α ^ 2 < 28 * d ^ 2 * Real.sqrt b := by
  have hX : 0 < 12 * d ^ 2 * b := by positivity
  have h2 : α ^ 2 ≤ 4 * Real.sqrt (12 * d ^ 2 * b) := by
    have : α ^ 2 ≤ (2 * (12 * d ^ 2 * b) ^ ((1 : ℝ) / 4)) ^ 2 := by gcongr
    rw [mul_pow, ← Real.rpow_natCast ((12 * d ^ 2 * b) ^ ((1 : ℝ) / 4)) 2,
      ← Real.rpow_mul hX.le, Real.sqrt_eq_rpow] at *
    norm_num at this ⊢; linarith
  have h3 : Real.sqrt (12 * d ^ 2 * b) = d * Real.sqrt 12 * Real.sqrt b := by
    rw [show 12 * d ^ 2 * b = d ^ 2 * (12 * b) by ring, Real.sqrt_mul (by positivity),
      Real.sqrt_sq hd.le, Real.sqrt_mul (by norm_num)]; ring
  have h12 : Real.sqrt 12 < 3.5 := by
    rw [Real.sqrt_lt' (by norm_num)]; norm_num
  have hsb : 0 < Real.sqrt b := Real.sqrt_pos.2 hb
  rw [h3] at h2
  have : 2 * d * α ^ 2 ≤ 8 * d ^ 2 * Real.sqrt 12 * Real.sqrt b := by nlinarith
  have : 8 * d ^ 2 * Real.sqrt 12 * Real.sqrt b < 28 * d ^ 2 * Real.sqrt b := by
    have : 0 < d ^ 2 * Real.sqrt b := by positivity
    nlinarith
  linarith

/-! ### §9.2.3 radix-digit polynomials -/

end IntMul.OpenAI.RP

end

section
-- Source module: Solutions.ResampleCost
/-!
# §9.8 The complete time bound

* `Ncost`: the normalized cost of Table 1 (cost divided by `V = T p`),
  `dK + pK^{τ-1} + ℓ d^{λ'} + d²(1+ℓ^τ) + d p^{1/2+δ} α + d p^δ + log₂(rp) + p^δ + 1`.
* `Ncost_le`: under the pointwise parameter bounds of §9.2.1
  (`p^ε/12 ≤ d ≤ p^ε`, `p^{εc}/24 ≤ K ≤ p^{εc}`, `ℓ ≤ p/(3d)`, `α ≤ 4 p^{1/4+ε/2}`),
  `Ncost ≤ 50 p^{1-2κ}` with `κ = 2⁻¹⁸²`.
* `concrete_d_bounds`, `concrete_K_bounds`, `concrete_alpha_bound`: the floors/ceilings
  `d = ⌊b^ε⌋`, `K = ⌊d^c⌋`, `α = ⌈(12d²b)^{1/4}⌉` satisfy these bounds with `p = 6b`.
* `kappaBound_of_table`: a correct multitape machine whose time on length `n ≥ n₀` is at most
  `C₀ (V(n) · Ncost(n) + n)` with `V ≤ 48n`, `p = 6 lg n` and parameters obeying the pointwise
  bounds satisfies the campaign statement `KappaBound (2⁻¹⁸²)`, i.e. time `O(n (lg n)^{1-κ})`.
-/

open Real

namespace IntMul.OpenAI.RCost

noncomputable section

/-! ### The concrete parameters of §9.2.1 satisfy the pointwise bounds -/

lemma floor_half {x : ℝ} (hx : 2 ≤ x) : x / 2 ≤ (⌊x⌋₊ : ℝ) := by
  have := Nat.lt_floor_add_one x
  linarith

/-- `d = ⌊b^ε⌋` with `b^ε ≥ 2` satisfies `p^ε/12 ≤ d ≤ p^ε` for `p = 6b`. -/
lemma concrete_d_bounds (b : ℝ) (hb : 1 ≤ b) (h2 : 2 ≤ b ^ ε) :
    (6 * b) ^ ε / 12 ≤ (⌊b ^ ε⌋₊ : ℝ) ∧ (⌊b ^ ε⌋₊ : ℝ) ≤ (6 * b) ^ ε := by
  have hb0 : 0 ≤ b := by linarith
  have h6 : (6 * b) ^ ε = 6 ^ ε * b ^ ε := Real.mul_rpow (by norm_num) hb0
  have h6le : (6 : ℝ) ^ ε ≤ 6 := by
    calc (6 : ℝ) ^ ε ≤ 6 ^ (1 : ℝ) :=
          Real.rpow_le_rpow_of_exponent_le (by norm_num) (by unfold ε; norm_num)
      _ = 6 := Real.rpow_one 6
  have h1le : (1 : ℝ) ≤ 6 ^ ε := Real.one_le_rpow (by norm_num) (by unfold ε; norm_num)
  have hbe : 0 ≤ b ^ ε := Real.rpow_nonneg hb0 ε
  constructor
  · rw [h6]
    have := floor_half h2
    nlinarith
  · rw [h6]
    have := Nat.floor_le hbe
    nlinarith

/-- `K = ⌊d^c⌋` with `d^c ≥ 2` and `p^ε/12 ≤ d ≤ p^ε` satisfies
`p^{εc}/24 ≤ K ≤ p^{εc}`. -/
lemma concrete_K_bounds (p d : ℝ) (hp : 1 ≤ p) (hd1 : p ^ ε / 12 ≤ d) (hd2 : d ≤ p ^ ε)
    (h2 : 2 ≤ d ^ c) :
    p ^ (ε * c) / 24 ≤ (⌊d ^ c⌋₊ : ℝ) ∧ (⌊d ^ c⌋₊ : ℝ) ≤ p ^ (ε * c) := by
  have hp0 : 0 < p := by linarith
  have hpe : 0 < p ^ ε := Real.rpow_pos_of_pos hp0 ε
  have hc0 : 0 ≤ c := by unfold c; norm_num
  have hd0 : 0 ≤ d := le_trans (by positivity) hd1
  constructor
  · have hfl := floor_half h2
    have hlow : (p ^ ε / 12) ^ c ≤ d ^ c := Real.rpow_le_rpow (by positivity) hd1 hc0
    have e : (p ^ ε / 12) ^ c = p ^ (ε * c) / 12 ^ c := by
      rw [Real.div_rpow hpe.le (by norm_num), ← Real.rpow_mul hp0.le]
    have h12 : (12 : ℝ) ^ c ≤ 12 := by
      calc (12 : ℝ) ^ c ≤ 12 ^ (1 : ℝ) :=
            Real.rpow_le_rpow_of_exponent_le (by norm_num) (by unfold c; norm_num)
        _ = 12 := Real.rpow_one 12
    have h12' : 0 < (12 : ℝ) ^ c := Real.rpow_pos_of_pos (by norm_num) c
    have : p ^ (ε * c) / 12 ≤ p ^ (ε * c) / 12 ^ c := by
      apply div_le_div_of_nonneg_left (Real.rpow_pos_of_pos hp0 _).le h12' h12
    rw [e] at hlow
    linarith
  · have := Nat.floor_le (Real.rpow_nonneg hd0 c)
    have hup : d ^ c ≤ (p ^ ε) ^ c := Real.rpow_le_rpow hd0 hd2 hc0
    rw [← Real.rpow_mul hp0.le] at hup
    linarith

/-- `α = ⌈(12d²b)^{1/4}⌉ ≤ 4 p^{1/4+ε/2}` for `p = 6b`, `d ≤ p^ε`, `(12d²b)^{1/4} ≥ 1`. -/
lemma concrete_alpha_bound (b d : ℝ) (hb : 1 ≤ b) (hd0 : 0 ≤ d) (hd : d ≤ (6 * b) ^ ε)
    (h1 : 1 ≤ (12 * d ^ 2 * b) ^ ((1 : ℝ) / 4)) :
    (⌈(12 * d ^ 2 * b) ^ ((1 : ℝ) / 4)⌉₊ : ℝ) ≤ 4 * (6 * b) ^ (1 / 4 + ε / 2) := by
  set X := (12 * d ^ 2 * b) ^ ((1 : ℝ) / 4)
  have hceil : (⌈X⌉₊ : ℝ) ≤ 2 * X := by
    have := Nat.ceil_lt_add_one (show 0 ≤ X by linarith)
    linarith
  have hp0 : 0 < 6 * b := by linarith
  have hpe : 0 < (6 * b) ^ ε := Real.rpow_pos_of_pos hp0 ε
  -- 12 d² b ≤ 2 (6b)^{2ε} (6b)
  have hY : 12 * d ^ 2 * b ≤ 2 * ((6 * b) ^ (2 * ε) * (6 * b)) := by
    have : d ^ 2 ≤ ((6 * b) ^ ε) ^ 2 := by gcongr
    have e : ((6 * b) ^ ε) ^ 2 = (6 * b) ^ (2 * ε) := by
      rw [← Real.rpow_natCast, ← Real.rpow_mul hp0.le]; norm_num; ring_nf
    rw [e] at this
    nlinarith
  have hX : X ≤ (2 * ((6 * b) ^ (2 * ε) * (6 * b))) ^ ((1 : ℝ) / 4) :=
    Real.rpow_le_rpow (by positivity) hY (by norm_num)
  have e2 : (2 * ((6 * b) ^ (2 * ε) * (6 * b))) ^ ((1 : ℝ) / 4) =
      2 ^ ((1 : ℝ) / 4) * (6 * b) ^ (1 / 4 + ε / 2) := by
    rw [Real.mul_rpow (by norm_num) (by positivity)]
    congr 1
    rw [← Real.rpow_add_one hp0.ne', ← Real.rpow_mul hp0.le]
    congr 1; ring
  have h2 : (2 : ℝ) ^ ((1 : ℝ) / 4) ≤ 2 := by
    calc (2 : ℝ) ^ ((1 : ℝ) / 4) ≤ 2 ^ (1 : ℝ) :=
          Real.rpow_le_rpow_of_exponent_le (by norm_num) (by norm_num)
      _ = 2 := Real.rpow_one 2
  have hq : 0 < (6 * b) ^ (1 / 4 + ε / 2) := Real.rpow_pos_of_pos hp0 _
  rw [e2] at hX
  nlinarith

/-! ### Connection to the campaign statement -/

end

end IntMul.OpenAI.RCost

end