Round-six κ development, definitions part 61

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

import Definitions.Def_IntMul_Kappa_R6_Defs60
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_RealExp
import Definitions.Def_IntMul_HvdH_Twiddle
import Definitions.Def_IntMul_HvdH_SApprox
import Definitions.Def_IntMul_HvdH_DiagApprox
import Definitions.Def_IntMul_HvdH_EApprox
import Definitions.Def_IntMul_HvdH_JApprox
import Definitions.Def_IntMul_HvdH_BApprox
import Definitions.Def_IntMul_HvdH_FPApprox
import Definitions.Def_IntMul_HvdH_MachineRef
import Definitions.Def_IntMul_Kappa_K_Defs9
import Definitions.Def_IntMul_Kappa_K_Defs19
import Definitions.Def_IntMul_Kappa_K_Defs14
import Definitions.Def_IntMul_Kappa_R6_Defs3

set_option autoImplicit false
set_option relaxedAutoImplicit false
set_option maxHeartbeats 1000000

section
-- Source module: Solutions.GlueParams
/-!
# Glue 3: explicit parameters `n ↦ (b, L, T, d, K, ℓ, α, p)` of §9.2 and the threshold `n₀`

For an input length `n` (§9.2.1):
* `b = lg n` (`bP`), `p = 6b` (`pP`), `d = ⌊b^ε⌋` (`dP`), `K = ⌊d^c⌋` (`KP`),
  `α = ⌈(12 d² b)^{1/4}⌉` (`αP`);
* `L` (`LP`) = the least `L` with `4n ≤ 2^L b`, `T = 2^L` (`TP`), so `4n ≤ T b < 8n`;
* `ℓ = ⌈L/d⌉` (`ℓP`, `RP.ell`), `r = 2^ℓ` (`rP`), bit volume `V = T p` (`VP`), `γ = 2dα²` (`γP`).

The threshold is `n₀ = 2^{b₀}` with `b₀ = 2^{2^{132}}` (`n0`, `b0`); for `n ≥ n₀` we have
`b ≥ b₀` (`b_ge`), and `params_ok` proves every pointwise condition used in §8–9:
`b^ε ≥ 2`, `d^c ≥ 2`, `1 ≤ d ≤ b`, `b/2 ≤ L < b`, `24n ≤ Tp ≤ 48n`, `p/(12d) ≤ ℓ ≤ p/(3d)`,
`64d²(256d²)! ≤ r`, `2 | r`, `2dp² < TL`, `b ≥ 8192`, `T < 2^b`, `12 d² b ≤ α⁴`,
`γ ≤ b/4`, the bounds of `kappaBound_of_table` on `d, K, ℓ, α`, and the existence of the
pairwise coprime odd source lengths (`RP.sources_exist`).

**`kappaBound_of_machine`**: the only hypotheses left are about the machine `M`
(correct for all `n ≥ 1`; time `≤ C₀ (V · Ncost + n)` for `n ≥ n₀` with these parameters).
-/

open Real

namespace IntMul.OpenAI.Glue

open RCost

noncomputable section

/-! ## 1. The parameter map -/

/-! ## 2. Elementary inequalities -/

lemma two_mul_succ_lt_pow : ∀ k, 3 ≤ k → 2 * k + 1 < 2 ^ k := by
  intro k hk
  induction k, hk using Nat.le_induction with
  | base => norm_num
  | succ k hk ih => rw [pow_succ]; omega

lemma cube_lt_pow : ∀ b, 40 ≤ b → 72 * b ^ 3 < 2 ^ b := by
  intro b hb
  induction b, hb using Nat.le_induction with
  | base => norm_num
  | succ b hb ih =>
    rw [pow_succ 2 b]
    have h1 : 40 * b ^ 2 ≤ b * b ^ 2 := Nat.mul_le_mul_right _ hb
    have h2 : 40 * b ≤ b * b := Nat.mul_le_mul_right _ hb
    have : (b + 1) ^ 3 ≤ 2 * b ^ 3 := by ring_nf; ring_nf at h1 h2; nlinarith
    nlinarith

lemma lin_lt_pow_half (b : ℕ) (hb : 16 ≤ b) : b < 2 ^ (b / 2) := by
  have := two_mul_succ_lt_pow (b / 2) (by omega)
  omega

/-- `64 x² (256 x²)! ≤ 2^{2^{17} x⁴}`. -/
lemma fact_bound (x : ℕ) : 64 * x ^ 2 * (256 * x ^ 2).factorial ≤ 2 ^ (2 ^ 17 * x ^ 4) := by
  set m := 256 * x ^ 2 with hm
  have h1 : m.factorial ≤ m ^ m := Nat.factorial_le_pow m
  have h2 : m ^ m ≤ (2 ^ m) ^ m := Nat.pow_le_pow_left (Nat.lt_two_pow_self).le m
  have h3 : 64 * x ^ 2 ≤ 2 ^ (64 * x ^ 2) := (Nat.lt_two_pow_self).le
  have h4 : 64 * x ^ 2 * m.factorial ≤ 2 ^ (64 * x ^ 2) * (2 ^ m) ^ m :=
    Nat.mul_le_mul h3 (h1.trans h2)
  rw [← pow_mul, ← pow_add] at h4
  refine h4.trans (Nat.pow_le_pow_right (by norm_num) ?_)
  rw [hm]
  nlinarith [sq_nonneg (x ^ 2), Nat.zero_le (x ^ 2), Nat.zero_le (x ^ 4)]

/-! ## 3. Consequences of `n ≥ n₀` -/

lemma b0_real_ge : (2 : ℝ) ^ ((2 : ℝ) ^ 132) ≤ (b0 : ℝ) := Nat.le_ceil _

lemma b0_ge (k : ℕ) (hk : k ≤ 2 ^ 132) : 2 ^ k ≤ b0 := by
  have hk' : (k : ℝ) ≤ (2 : ℝ) ^ 132 := by exact_mod_cast hk
  have : (2 : ℝ) ^ k ≤ (b0 : ℝ) := by
    rw [← Real.rpow_natCast]
    exact (Real.rpow_le_rpow_of_exponent_le (by norm_num) hk').trans b0_real_ge
  exact_mod_cast this

lemma b_ge {n : ℕ} (hn : n0 ≤ n) : b0 ≤ bP n := by
  have h := Nat.clog_mono_right 2 hn
  rw [n0, Nat.clog_pow 2 b0 (by norm_num)] at h
  exact h.trans (le_max_left _ _)

lemma bP_eq {n : ℕ} (hn : n0 ≤ n) : bP n = Nat.clog 2 n := by
  have h := Nat.clog_mono_right 2 hn
  rw [n0, Nat.clog_pow 2 b0 (by norm_num)] at h
  have : 1 ≤ b0 := b0_ge 0 (by norm_num)
  exact max_eq_left (by omega)

lemma n_le {n : ℕ} (hn : n0 ≤ n) : n ≤ 2 ^ bP n := by
  rw [bP_eq hn]; exact Nat.le_pow_clog (by norm_num) n

lemma n_gt {n : ℕ} (hn : n0 ≤ n) : 2 ^ (bP n - 1) < n := by
  have h2 : 2 ≤ n0 := by
    have : 2 ^ 1 ≤ 2 ^ b0 := Nat.pow_le_pow_right (by norm_num) (b0_ge 0 (by norm_num))
    simpa [n0] using this
  rw [bP_eq hn]
  have := Nat.pow_pred_clog_lt_self (b := 2) (by norm_num) (show 1 < n by omega)
  simpa [Nat.pred_eq_sub_one] using this

/-- Facts about `L` and `T` (`b ≥ 16`). -/
lemma L_facts {n : ℕ} (hn : n0 ≤ n) :
    4 * n ≤ TP n * bP n ∧ TP n * bP n < 8 * n ∧ LP n < bP n ∧ bP n ≤ 2 * LP n := by
  have hb : 2 ^ 4 ≤ bP n := (b0_ge 4 (by norm_num)).trans (b_ge hn)
  have hb' : 16 ≤ bP n := by simpa using hb
  have hle := n_le hn
  have hgt := n_gt hn
  set b := bP n with hbdef
  have hpow : 2 ^ b = 2 * 2 ^ (b - 1) := by
    rw [← pow_succ']; congr 1; omega
  have spec : 4 * n ≤ 2 ^ LP n * b := Nat.find_spec (LP_exists n)
  -- `L < b`
  have hLb : LP n ≤ b - 1 := by
    apply Nat.find_min' (LP_exists n)
    show 4 * n ≤ 2 ^ (b - 1) * b
    have : 4 * n ≤ 8 * 2 ^ (b - 1) := by omega
    nlinarith
  -- `b ≤ 2L`
  have h2L : b ≤ 2 * LP n := by
    by_contra hcon
    push Not at hcon
    have hhalf : LP n + b / 2 ≤ b := by omega
    have h1 : b < 2 ^ (b / 2) := lin_lt_pow_half b hb'
    have h3 : 2 ^ LP n * b < 2 ^ LP n * 2 ^ (b / 2) := Nat.mul_lt_mul_of_pos_left h1 (by positivity)
    have h4 : 2 ^ LP n * 2 ^ (b / 2) ≤ 2 ^ b := by
      rw [← pow_add]; exact Nat.pow_le_pow_right (by norm_num) hhalf
    omega
  -- `T b < 8n`
  have hTb : TP n * b < 8 * n := by
    have hL1 : 1 ≤ LP n := by omega
    have hmin : ¬ 4 * n ≤ 2 ^ (LP n - 1) * b := Nat.find_min (LP_exists n) (m := LP n - 1)
      (by show LP n - 1 < LP n; omega)
    have : TP n = 2 * 2 ^ (LP n - 1) := by
      rw [TP, ← pow_succ']; congr 1; omega
    rw [this]; push Not at hmin; nlinarith
  exact ⟨spec, hTb, by omega, h2L⟩

lemma beps_big {n : ℕ} (hn : n0 ≤ n) : (2 : ℝ) ^ ((2 : ℝ) ^ 57) ≤ (bP n : ℝ) ^ ε := by
  have hB : (2 : ℝ) ^ ((2 : ℝ) ^ 132) ≤ (bP n : ℝ) := by
    exact b0_real_ge.trans (by exact_mod_cast b_ge hn)
  have h0 : 0 ≤ (2 : ℝ) ^ ((2 : ℝ) ^ 132) := by positivity
  calc (2 : ℝ) ^ ((2 : ℝ) ^ 57) = ((2 : ℝ) ^ ((2 : ℝ) ^ 132)) ^ ε := by
        rw [← Real.rpow_mul (by norm_num), show (2 : ℝ) ^ 132 * ε = (2 : ℝ) ^ 57 by
          unfold ε; norm_num]
    _ ≤ (bP n : ℝ) ^ ε := Real.rpow_le_rpow h0 hB (by unfold ε; norm_num)

lemma two_le_beps {n : ℕ} (hn : n0 ≤ n) : 2 ≤ (bP n : ℝ) ^ ε := by
  refine le_trans ?_ (beps_big hn)
  calc (2 : ℝ) = 2 ^ (1 : ℝ) := (Real.rpow_one 2).symm
    _ ≤ (2 : ℝ) ^ ((2 : ℝ) ^ 57) := Real.rpow_le_rpow_of_exponent_le (by norm_num) (by norm_num)

lemma d_half {n : ℕ} (hn : n0 ≤ n) : (bP n : ℝ) ^ ε / 2 ≤ (dP n : ℝ) :=
  RCost.floor_half (two_le_beps hn)

lemma d_le_beps (n : ℕ) : (dP n : ℝ) ≤ (bP n : ℝ) ^ ε :=
  Nat.floor_le (Real.rpow_nonneg (Nat.cast_nonneg _) _)

lemma two_le_dc {n : ℕ} (hn : n0 ≤ n) : 2 ≤ (dP n : ℝ) ^ c := by
  have hY : (2 : ℝ) ^ ((2 : ℝ) ^ 57 - 1) ≤ (dP n : ℝ) := by
    rw [Real.rpow_sub_one (by norm_num)]
    exact le_trans (div_le_div_of_nonneg_right (beps_big hn) (by norm_num)) (d_half hn)
  have h0 : 0 ≤ (2 : ℝ) ^ ((2 : ℝ) ^ 57 - 1) := by positivity
  calc (2 : ℝ) = 2 ^ (1 : ℝ) := (Real.rpow_one 2).symm
    _ ≤ (2 : ℝ) ^ (((2 : ℝ) ^ 57 - 1) * c) :=
        Real.rpow_le_rpow_of_exponent_le (by norm_num) (by unfold c; norm_num)
    _ = ((2 : ℝ) ^ ((2 : ℝ) ^ 57 - 1)) ^ c := by rw [Real.rpow_mul (by norm_num)]
    _ ≤ (dP n : ℝ) ^ c := Real.rpow_le_rpow h0 hY (by unfold c; norm_num)

lemma d_pos {n : ℕ} (hn : n0 ≤ n) : 1 ≤ dP n := by
  have h1 := d_half hn
  have h2 := two_le_beps hn
  have : (1 : ℝ) ≤ dP n := by linarith
  exact_mod_cast this

/-- `B ≥ 2^{36}`-type bounds from `b ≥ b₀`. -/
lemma b_ge_real {n : ℕ} (hn : n0 ≤ n) (k : ℕ) (hk : k ≤ 2 ^ 132) : (2 : ℝ) ^ k ≤ (bP n : ℝ) := by
  have := (b0_ge k hk).trans (b_ge hn)
  exact_mod_cast this

lemma d_le_b {n : ℕ} (hn : n0 ≤ n) : dP n ≤ bP n := by
  have hb1 : (1 : ℝ) ≤ bP n := by exact_mod_cast bP_pos n
  have : (dP n : ℝ) ≤ bP n := (d_le_beps n).trans (by
    calc (bP n : ℝ) ^ ε ≤ (bP n : ℝ) ^ (1 : ℝ) :=
          Real.rpow_le_rpow_of_exponent_le hb1 (by unfold ε; norm_num)
      _ = bP n := Real.rpow_one _)
  exact_mod_cast this

/-- `x^e ≥ 2^{k e}` style: `(2^k)^{1/2} = 2^{k/2}`. -/
lemma sqrt_pow_two (k : ℕ) : ((2 : ℝ) ^ (2 * k)) ^ ((1 : ℝ) / 2) = 2 ^ k := by
  rw [← Real.rpow_natCast, ← Real.rpow_mul (by norm_num), ← Real.rpow_natCast]
  congr 1; push_cast; ring

/-- `b ≥ 2^{18} d^5`. -/
lemma b_ge_d5 {n : ℕ} (hn : n0 ≤ n) : (2 : ℝ) ^ 18 * (dP n : ℝ) ^ 5 ≤ bP n := by
  set B : ℝ := (bP n : ℝ)
  have hB : (2 : ℝ) ^ 36 ≤ B := b_ge_real hn 36 (by norm_num)
  have hB1 : 1 ≤ B := le_trans (by norm_num) hB
  have hB0 : 0 ≤ B := by linarith
  have hd5 : (dP n : ℝ) ^ 5 ≤ B ^ (5 * ε) := by
    calc (dP n : ℝ) ^ 5 ≤ (B ^ ε) ^ 5 := by
          gcongr; exact d_le_beps n
      _ = B ^ (5 * ε) := by
          rw [← Real.rpow_natCast, ← Real.rpow_mul hB0]; congr 1; push_cast; ring
  have hsq : (2 : ℝ) ^ 18 ≤ B ^ ((1 : ℝ) / 2) := by
    rw [← sqrt_pow_two 18]
    exact Real.rpow_le_rpow (by positivity) (by simpa using hB) (by norm_num)
  have hhalf : B ^ ((1 : ℝ) / 2) ≤ B ^ (1 - 5 * ε) :=
    Real.rpow_le_rpow_of_exponent_le hB1 (by unfold ε; norm_num)
  have hsplit : B = B ^ (1 - 5 * ε) * B ^ (5 * ε) := by
    rw [← Real.rpow_add (by linarith)]; simp
  have hpos : 0 ≤ (dP n : ℝ) ^ 5 := by positivity
  calc (2 : ℝ) ^ 18 * (dP n : ℝ) ^ 5 ≤ B ^ (1 - 5 * ε) * B ^ (5 * ε) := by
        apply mul_le_mul (hsq.trans hhalf) hd5 hpos (Real.rpow_nonneg hB0 _)
    _ = B := hsplit.symm

/-- `28 d² √b ≤ b/4`. -/
lemma gamma_room {n : ℕ} (hn : n0 ≤ n) : 28 * (dP n : ℝ) ^ 2 * Real.sqrt (bP n) ≤ (bP n : ℝ) / 4 := by
  set B : ℝ := (bP n : ℝ)
  have hB : (2 : ℝ) ^ 36 ≤ B := b_ge_real hn 36 (by norm_num)
  have hB1 : 1 ≤ B := le_trans (by norm_num) hB
  have hB0 : 0 ≤ B := by linarith
  have hd2 : (dP n : ℝ) ^ 2 ≤ B ^ (2 * ε) := by
    calc (dP n : ℝ) ^ 2 ≤ (B ^ ε) ^ 2 := by gcongr; exact d_le_beps n
      _ = B ^ (2 * ε) := by
          rw [← Real.rpow_natCast, ← Real.rpow_mul hB0]; congr 1; push_cast; ring
  -- `112 B^{2ε} √B ≤ B`
  have hq : (2 : ℝ) ^ 9 ≤ B ^ ((1 : ℝ) / 4) := by
    have : ((2 : ℝ) ^ 36) ^ ((1 : ℝ) / 4) = 2 ^ 9 := by
      rw [← Real.rpow_natCast, ← Real.rpow_mul (by norm_num), ← Real.rpow_natCast]
      congr 1; norm_num
    rw [← this]; exact Real.rpow_le_rpow (by positivity) hB (by norm_num)
  have hq2 : B ^ ((1 : ℝ) / 4) ≤ B ^ (1 / 2 - 2 * ε) :=
    Real.rpow_le_rpow_of_exponent_le hB1 (by unfold ε; norm_num)
  have hsplit : B = B ^ (2 * ε) * B ^ ((1 : ℝ) / 2) * B ^ (1 / 2 - 2 * ε) := by
    rw [← Real.rpow_add (by linarith), ← Real.rpow_add (by linarith)]
    rw [show 2 * ε + 1 / 2 + (1 / 2 - 2 * ε) = (1 : ℝ) by ring, Real.rpow_one]
  have hsqrt : Real.sqrt B = B ^ ((1 : ℝ) / 2) := Real.sqrt_eq_rpow B
  rw [hsqrt]
  have h1 : 0 ≤ B ^ ((1 : ℝ) / 2) := Real.rpow_nonneg hB0 _
  have h2 : 0 ≤ B ^ (2 * ε) := Real.rpow_nonneg hB0 _
  have h112 : (112 : ℝ) ≤ B ^ (1 / 2 - 2 * ε) := le_trans (by norm_num) (hq.trans hq2)
  calc 28 * (dP n : ℝ) ^ 2 * B ^ ((1 : ℝ) / 2) ≤ 28 * B ^ (2 * ε) * B ^ ((1 : ℝ) / 2) := by
        gcongr
    _ = (B ^ (2 * ε) * B ^ ((1 : ℝ) / 2) * 112) / 4 := by ring
    _ ≤ (B ^ (2 * ε) * B ^ ((1 : ℝ) / 2) * B ^ (1 / 2 - 2 * ε)) / 4 := by gcongr
    _ = B / 4 := by rw [← hsplit]

/-! ## 4. All parameter conditions for `n ≥ n₀` -/

/-- **The pointwise parameter conditions of §8–9 hold for `n ≥ n₀`.** -/
theorem params_ok {n : ℕ} (hn : n0 ≤ n) :
    -- thresholds
    2 ≤ (bP n : ℝ) ^ ε ∧ 2 ≤ (dP n : ℝ) ^ c ∧ 8192 ≤ bP n ∧ 1 ≤ dP n ∧ dP n ≤ bP n ∧
    -- box exponent and volume
    bP n ≤ 2 * LP n ∧ LP n < bP n ∧ TP n < 2 ^ bP n ∧
    24 * (n : ℝ) ≤ VP n ∧ VP n ≤ 48 * n ∧
    -- axis exponent, factorial bound, source lengths
    (pP n : ℝ) / (12 * dP n) ≤ ℓP n ∧ (ℓP n : ℝ) ≤ (pP n : ℝ) / (3 * dP n) ∧
    64 * dP n ^ 2 * (256 * dP n ^ 2).factorial ≤ rP n ∧ 2 ∣ rP n ∧
    -- resampling / error budget
    2 * dP n * pP n ^ 2 < TP n * LP n ∧ 12 * (dP n : ℝ) ^ 2 * bP n ≤ (αP n : ℝ) ^ 4 ∧
    γP n ≤ (bP n : ℝ) / 4 ∧
    -- the hypotheses of `kappaBound_of_table`
    (6 * (IntMul.lg n : ℝ)) ^ ε / 12 ≤ dP n ∧ (dP n : ℝ) ≤ (6 * (IntMul.lg n : ℝ)) ^ ε ∧
    (6 * (IntMul.lg n : ℝ)) ^ (ε * c) / 24 ≤ KP n ∧ (KP n : ℝ) ≤ (6 * (IntMul.lg n : ℝ)) ^ (ε * c) ∧
    (ℓP n : ℝ) ≤ 6 * (IntMul.lg n : ℝ) / (3 * dP n) ∧
    (αP n : ℝ) ≤ 4 * (6 * (IntMul.lg n : ℝ)) ^ (1 / 4 + ε / 2) := by
  obtain ⟨hT1, hT2, hLb, h2L⟩ := L_facts hn
  have hbe := two_le_beps hn
  have hdc := two_le_dc hn
  have hd1 := d_pos hn
  have hdb := d_le_b hn
  have hb8 : 8192 ≤ bP n := le_trans (by norm_num) ((b0_ge 13 (by norm_num)).trans (b_ge hn))
  have hb1 : (1 : ℝ) ≤ bP n := by exact_mod_cast bP_pos n
  -- d and K bounds
  obtain ⟨hD1, hD2⟩ := concrete_d_bounds (bP n : ℝ) hb1 hbe
  have hp1 : (1 : ℝ) ≤ 6 * bP n := by linarith
  obtain ⟨hK1, hK2⟩ := concrete_K_bounds (6 * bP n) (dP n) hp1 hD1 hD2 hdc
  -- ℓ bounds
  obtain ⟨hl1, hl2⟩ := RP.ell_bounds (LP n) (dP n) (bP n) (by omega) hdb h2L hLb.le
  -- α
  have hX1 : (1 : ℝ) ≤ (12 * (dP n : ℝ) ^ 2 * bP n) ^ ((1 : ℝ) / 4) := by
    apply Real.one_le_rpow _ (by norm_num)
    have : (1 : ℝ) ≤ dP n := by exact_mod_cast hd1
    nlinarith
  have hα := concrete_alpha_bound (bP n) (dP n) hb1 (Nat.cast_nonneg _) hD2 hX1
  have hα4 : 12 * (dP n : ℝ) ^ 2 * bP n ≤ (αP n : ℝ) ^ 4 := by
    have h0 : 0 ≤ 12 * (dP n : ℝ) ^ 2 * bP n := by positivity
    have hc : (12 * (dP n : ℝ) ^ 2 * bP n) ^ ((1 : ℝ) / 4) ≤ αP n := Nat.le_ceil _
    have : ((12 * (dP n : ℝ) ^ 2 * bP n) ^ ((1 : ℝ) / 4)) ^ 4 = 12 * (dP n : ℝ) ^ 2 * bP n := by
      rw [← Real.rpow_natCast, ← Real.rpow_mul h0]; norm_num
    rw [← this]; gcongr
  -- γ
  have hγ : γP n ≤ (bP n : ℝ) / 4 := by
    have hαle : (αP n : ℝ) ≤ 2 * (12 * (dP n : ℝ) ^ 2 * bP n) ^ ((1 : ℝ) / 4) := by
      have := Nat.ceil_lt_add_one (show 0 ≤ (12 * (dP n : ℝ) ^ 2 * bP n) ^ ((1 : ℝ) / 4) by
        positivity)
      unfold αP; linarith
    have hd0 : (0 : ℝ) < dP n := by exact_mod_cast hd1
    have := RP.gamma_bound (dP n) (bP n) (αP n) hd0 (by linarith) (Nat.cast_nonneg _) hαle
    unfold γP; linarith [gamma_room hn]
  -- factorial bound
  have hfact : 64 * dP n ^ 2 * (256 * dP n ^ 2).factorial ≤ rP n := by
    refine (fact_bound (dP n)).trans (Nat.pow_le_pow_right (by norm_num) ?_)
    have h5 := b_ge_d5 hn
    have hd0 : (0 : ℝ) < dP n := by exact_mod_cast hd1
    have : ((2 ^ 17 * dP n ^ 4 : ℕ) : ℝ) ≤ ℓP n := by
      refine le_trans ?_ hl1
      push_cast
      rw [le_div_iff₀ (by positivity)]
      nlinarith
    exact_mod_cast this
  have hr2 : 2 ∣ rP n := by
    have hℓ1 : 1 ≤ ℓP n := by
      have : (1 : ℝ) ≤ ℓP n := by
        refine le_trans ?_ hl1
        rw [le_div_iff₀ (by positivity)]
        have h5 := b_ge_d5 hn
        have hd' : (1 : ℝ) ≤ dP n := by exact_mod_cast hd1
        have : (dP n : ℝ) ≤ (dP n : ℝ) ^ 5 := le_self_pow₀ hd' (by norm_num)
        nlinarith
      exact_mod_cast this
    exact dvd_pow_self 2 (by omega)
  -- 2 d p² < T L
  have hdp : 2 * dP n * pP n ^ 2 < TP n * LP n := by
    have hcube := cube_lt_pow (bP n) (by omega)
    have h1 : 2 * dP n * pP n ^ 2 ≤ 72 * bP n ^ 3 := by
      unfold pP; nlinarith
    have hgt := n_gt hn
    have hpow : 2 ^ bP n = 2 * 2 ^ (bP n - 1) := by rw [← pow_succ']; congr 1; omega
    have h2 : 2 ^ bP n < TP n * LP n * 2 := by nlinarith
    have h3 : 2 ^ bP n ≤ TP n * LP n := by
      have : TP n * bP n ≤ TP n * (2 * LP n) := Nat.mul_le_mul_left _ h2L
      nlinarith
    omega
  -- volume
  have hV1 : 24 * (n : ℝ) ≤ VP n := by
    unfold VP pP; push_cast
    have : ((4 * n : ℕ) : ℝ) ≤ ((TP n * bP n : ℕ) : ℝ) := by exact_mod_cast hT1
    push_cast at this; nlinarith
  have hV2 : VP n ≤ 48 * n := by
    unfold VP pP; push_cast
    have : ((TP n * bP n : ℕ) : ℝ) ≤ ((8 * n : ℕ) : ℝ) := by exact_mod_cast hT2.le
    push_cast at this; nlinarith
  have hTlt : TP n < 2 ^ bP n := Nat.pow_lt_pow_right (by norm_num) hLb
  refine ⟨hbe, hdc, hb8, hd1, hdb, h2L, hLb, hTlt, hV1, hV2, ?_, ?_, hfact, hr2, hdp, hα4, hγ,
    hD1, hD2, ?_, ?_, ?_, hα⟩
  · unfold pP; push_cast; exact hl1
  · unfold pP; push_cast; exact hl2
  · exact hK1
  · exact hK2
  · exact hl2

/-! ## 5. The final reduction -/

end

end IntMul.OpenAI.Glue

end

section
-- Source module: Solutions.KFMErr
/-!
# KFM 14: the error chain of the computed pipeline (§9.7.1), abstract form

Arrays are functions `ι → ℂ` with the sup norm; `hp = 2^{-p}`. Every computed procedure is
an arbitrary function with an error bound on grid inputs (`Grid p f`: `Q_p f = f`), so the IL
denotations (and `P18IL`) plug in later.

* `QA p f`: `Q_p` on every entry; `QA_grid`, `norm_QA_le`, `norm_QA_sub_le` (`≤ √2 hp`).
* **`phase_err`**: `‖Q_p(ẽ·x) - e·x‖ ≤ 4hp`, `‖Q_p(ẽ·x)‖ ≤ ‖x‖` (unit phases, `|ẽ-e| ≤ 2hp`).
* `BilC M`: bilinear contraction; **`bil_err`**.
* **`mc_err`**: the computed complex convolution `Q_p(θ̃⁻ · M̃r(Q_p(θ̃⁺f), Q_p(θ̃⁺g)))` is within
  `(Er + 12) hp` of `θ⁻ · Mr(θ⁺ f, θ⁺ g)` (one operand of norm `≤ .7`).
* **`chirp_err`**: `Q_p(ã · M̃C(ã, Q_p(ã · x)))` is within `(Ec + 10) hp` of `a · MC(a, a · x)`.
* **`source_err`**: `2^γ B̃₀(C̃(Ã u))` is within `2^γ (εA + Ec + εB) hp` of `2^γ B₀(QF(A u))`.
* **`pipeline_err`**: with `z = Q_p(F̃u · F̃v)` and `F̃⁺ = conj ∘ F̃ ∘ conj`,
  `‖F̃⁺ z - G(Fu · Fv)‖ ≤ (3Es + 2) hp`.
-/

namespace NLogN.KFM

open Complex

noncomputable section

/-! ## Bilinear contractions -/

section Bil

variable {ι κ : Type*} [Fintype ι] [Fintype κ]

lemma LinC.err {L : (ι → ℂ) → (κ → ℂ)} (hL : LinC L) (a b : ι → ℂ) : ‖L a - L b‖ ≤ ‖a - b‖ := by
  rw [hL.sub]; exact hL.le _

end Bil

/-! ## The computed complex convolution `M̃_C` (§9.4.1) -/

/-! ## The computed chirp transform (§9.5) -/

/-! ## The source transform (§9.5.1) -/

/-! ## The pipeline (§9.7.2) -/

end

end NLogN.KFM

end

section
-- Source module: Solutions.KFMSetup
/-!
# KFM 14: the setup arithmetic (§9.2), in the integer form the machine computes

* **`floor_root_eq`/`floor_root_spec`**: `⌊N^{1/m}⌋ = x` iff `x^m ≤ N < (x+1)^m`, so
  **`dP_spec`** (`m = 2^75`, of `b`) and **`KP_spec`** (`m = 2^56`, of `d`) are exact integer
  comparisons (computed by saturating powers and binary search).
* **`ceil_root_eq`**, **`αP_spec`**: `α = αP n` iff `(α-1)^4 < 12 d² b ≤ α^4`.
* **`LP_spec`/`LP_eq`**: `L = LP n` iff `4n ≤ 2^L b` and (`L = 0` or `4n > 2^{L-1} b`).
* **`srcS d d' r i`**: the explicit source lengths of `RCop.coprime_params`
  (`1 + Q m_i`, `Q = 2F(⌊r/(2FD)⌋+1)`, `F = D!`, `D = 256 d²`), and **`srcS_ok`**: odd, pairwise
  coprime, in the windows `((1-2η)t_i, (1-η)t_i]`; **`srcS_ok_n`** for the parameters of `n ≥ n₀`.
* **`mu_spec`**: `μ = (P mod s)⁻¹` from `Nat.gcdA` (extended Euclid) satisfies `μ P ≡ 1 (mod s)`.
-/

namespace NLogN.KFM

open IntMul.OpenAI.Glue IntMul.OpenAI.RCost

noncomputable section

/-! ## Integer roots -/

lemma root_le_iff (N m x : ℕ) (hm : 1 ≤ m) :
    (x : ℝ) ≤ (N : ℝ) ^ ((1 : ℝ) / m) ↔ x ^ m ≤ N := by
  have hm' : (0 : ℝ) < m := by exact_mod_cast hm
  rw [one_div, Real.le_rpow_inv_iff_of_pos (Nat.cast_nonneg _) (Nat.cast_nonneg _) hm',
    Real.rpow_natCast]
  exact_mod_cast Iff.rfl

theorem floor_root_spec (N m : ℕ) (hm : 1 ≤ m) :
    ⌊(N : ℝ) ^ ((1 : ℝ) / m)⌋₊ ^ m ≤ N ∧ N < (⌊(N : ℝ) ^ ((1 : ℝ) / m)⌋₊ + 1) ^ m := by
  set x := ⌊(N : ℝ) ^ ((1 : ℝ) / m)⌋₊
  refine ⟨(root_le_iff N m x hm).1 (Nat.floor_le (by positivity)), ?_⟩
  by_contra h
  push Not at h
  have h3 := (root_le_iff N m (x + 1) hm).2 h
  have h4 := Nat.lt_floor_add_one ((N : ℝ) ^ ((1 : ℝ) / m))
  push_cast at h3
  linarith

lemma eps_eq : ε = (1 : ℝ) / ((2 ^ 75 : ℕ) : ℝ) := by unfold ε; norm_num

/-- `d = dP n` is the integer `2^75`-th root of `b`. -/
theorem dP_spec_r6c (n : ℕ) : dP n ^ (2 ^ 75) ≤ bP n ∧ bP n < (dP n + 1) ^ (2 ^ 75) := by
  have := floor_root_spec (bP n) (2 ^ 75) Nat.one_le_two_pow
  unfold dP; rw [eps_eq]; exact this

/-! ## Explicit source lengths -/

/-! ## Modular inverses -/

end

end NLogN.KFM

end

section
-- Source module: Solutions.KFMS4
/-!
# KFM S4: the cost of the scalar setup program

For `n ≥ n₀`, with `b = lg n`, every atom the setup expressions read is at most `A = 2^{b+2}`.
The expressions have degree at most `8`, so all values have fewer than `B = 8b+17` bits
(`KFMSm`). This gives

* **`setupA_cost`**, **`setupB_cost`**, **`setupC_cost`**: each at most `10^15 (b+1)^4`.
-/

set_option linter.unusedSimpArgs false
set_option linter.unusedVariables false

namespace NLogN.KFM

open IntMul IntMul.OpenAI.Glue IntMul.OpenAI.RP RS

/-- The size facts used by the cost bounds. -/
theorem setup_facts_r6c {n : ℕ} (hn : n0 ≤ n) :
    8192 ≤ bP n ∧ 1 ≤ dP n ∧ dP n ≤ bP n ∧ dP n ^ 2 ≤ bP n ∧ n ≤ 2 ^ bP n ∧ LP n < bP n ∧
    ℓP n ≤ LP n ∧ 64 * dP n ^ 2 * (256 * dP n ^ 2).factorial ≤ rP n ∧ rP n ≤ 2 ^ bP n ∧
    12 * dP n ^ 2 * bP n ≤ 2 ^ bP n ∧ 64 ≤ rP n ∧ srcQ (dP n) (rP n) ≤ 3 * 2 ^ bP n := by
  obtain ⟨-, -, hb, hd1, hdb, -, hLb, -, -, -, -, -, hbig, -⟩ := params_ok hn
  have hd2 : dP n ^ 2 ≤ bP n :=
    le_trans (Nat.pow_le_pow_right hd1 (by norm_num : 2 ≤ 2 ^ 75)) (dP_spec_r6c n).1
  have hnb : n ≤ 2 ^ bP n :=
    le_trans (Nat.le_pow_clog (by norm_num) n) (Nat.pow_le_pow_right (by norm_num) (le_max_left _ _))
  have hel : ℓP n ≤ LP n := by
    obtain ⟨-, h2⟩ := ell_spec (LP n) (dP n) (by omega)
    unfold ℓP
    by_contra hc
    push Not at hc
    have : dP n * (LP n + 1) ≤ dP n * ell (LP n) (dP n) := Nat.mul_le_mul_left _ hc
    rw [Nat.mul_add, mul_one] at this
    have : LP n ≤ dP n * LP n := Nat.le_mul_of_pos_left _ hd1
    omega
  have hr : rP n ≤ 2 ^ bP n := Nat.pow_le_pow_right (by norm_num) (by omega)
  have hF1 : 1 ≤ (256 * dP n ^ 2).factorial := Nat.one_le_iff_ne_zero.2 (Nat.factorial_ne_zero _)
  have h1d : 1 ≤ dP n ^ 2 := Nat.one_le_pow _ _ hd1
  have hr64 : 64 ≤ rP n := by nlinarith
  have hX : 12 * dP n ^ 2 * bP n ≤ 2 ^ bP n := by
    have := cube_lt_pow (bP n) (by omega)
    have : bP n ≤ bP n ^ 2 := by nlinarith
    nlinarith
  refine ⟨hb, hd1, hdb, hd2, hnb, hLb, hel, hbig, hr, hX, hr64, ?_⟩
  -- Q ≤ r + 2F
  unfold srcQ
  set F := (256 * dP n ^ 2).factorial
  set D := 256 * dP n ^ 2
  have hD1 : 1 ≤ D := by omega
  have h1 : 2 * F * (rP n / (2 * F * D)) ≤ rP n := by
    calc 2 * F * (rP n / (2 * F * D)) ≤ 2 * F * D * (rP n / (2 * F * D)) :=
          Nat.mul_le_mul_right _ (Nat.le_mul_of_pos_right _ hD1)
      _ ≤ rP n := Nat.mul_div_le _ _
  have h2 : 2 * F ≤ rP n := by nlinarith
  rw [Nat.mul_add, mul_one]; omega

end NLogN.KFM

end

section
-- Source module: Solutions.J6AMath4
/-!
# A lines (PR #5 block chirp), part 4: the generic `A` pass

`passGen` is `KFY1Pad.passA` with the per-axis fixed-point maps `Aa i` arbitrary: they must read
only the first `s_i` entries of a line (`hloc`) and be `FPApprox`es of `liftL (L i)` within `ε`
for line maps `‖L i‖ ≤ 1`. The list pipeline is `tensorGo (mapsGen …)`.
-/

set_option linter.unusedSimpArgs false
set_option linter.unusedVariables false
set_option linter.unusedSectionVars false

namespace J6A

open IntMul IntMul.HvdH IntMul.HvdH.Ref IntMul.HvdH.RefPf NLogN.KFM NLogN.KFY IntMul.OpenAI.RT
open Complex
open NLogN.KFB (GI gflat)

section gen
variable {d : ℕ} (p : ℕ) (s t : Fin d → ℕ) [∀ i, NeZero (s i)] [∀ i, NeZero (t i)]

/-- The `tensorGo` maps of per-axis line functions `Aa i` (axis `i`: `t_i` outputs). -/
def mapsGen (Aa : ℕ → (ℕ → ℤ × ℤ) → (ℕ → ℤ × ℤ)) : List (ℕ × (List GI → List GI)) :=
  (List.range d).map fun i => ((List.ofFn t).getD i 1, fun X : List (ℤ × ℤ) =>
    (List.range ((List.ofFn t).getD i 1)).map (Aa i fun j => X.getD j (0, 0)))

end gen

end J6A

end

section
-- Source module: Solutions.J6AMath5
/-!
# A lines (PR #5 block chirp), part 5: the configured chirp map

Precision `P = 2p`, output exponent `σ = ⌊p/2⌋`, half-block `1 ≤ H ≤ 2m` (`m = sWidth p α`),
window half-width `W = H + m + 3`.

* `chirpA_congr`: `chirpA` reads only `[0, s)`.
* `chirp_sigma`: `exp(π s(t-s)H²/(t²α²)) ≤ 2^{⌊p/2⌋}` once `300 (t-s) ≤ s`.
* **`chirpC_fp`**: `chirpC` is an `FPApprox` of `liftL L` within `12 m + 11` for every line map
  `L` fitted to `Sapprox` within `6m+4` (the `FitSpec6` clause).
-/

set_option linter.unusedVariables false
set_option linter.unusedSectionVars false

open Complex

namespace J6A

open IntMul.HvdH

/-- The configured chirp line map. -/
def chirpC (p s t α H : ℕ) : (ℕ → ℤ × ℤ) → (ℕ → ℤ × ℤ) :=
  chirpA (2 * p) (p / 2) (H + sWidth p α + 3) s t α H

end J6A

end

section
-- Source module: Solutions.J6SegGlob
/-!
# Segmented inverse (S2): the resampling system split into blocks, coupling and tail

For a valid block structure (`BlocksOK`: blocks are wrap-free runs `[st i, st i + len i)` that
partition `[0, s)`, coupling width `w ≤ s`):
* `MxG` (the block part of `N`, entries `δ_{ik} + e_{i,k-i}`) and `MiG` (its exact block inverse
  `D T⁻¹ D⁻¹` from the q-LU) satisfy `M M⁻¹ = I` (`hinvG`, via `eEx_recentre`, `MF_mul_MinvF`);
* row sums: `‖M - I‖, ‖Δ‖ ≤ 2^{-α²θ}` (`rowM`, `rowD`), hence `‖M⁻¹‖ ≤ 8/7` when `α²θ ≥ 3`;
* every coupling column lies in the capacitance set `Jlist` (`col_mem`);
* `nsplit`: `(N x)_i = (M x)_i + (Δ x)_i + drop`, `|drop| ≤ 4 e^{-πα² w(w+1)} ‖x‖`;
* `nearN`: `‖N d‖ ≤ ρ ⇒ ‖d‖ ≤ 8ρ/7` (`‖N - I‖ ≤ 1/8`).
-/

open Finset Real

namespace J6Seg
open IntMul.HvdH IntMul.HvdH.JPf J6QLU

section glob
variable {s t α : ℕ} {B : Blocks}

lemma toNat_emod_lt (hs : 0 < s) (a : ℤ) : (a % s).toNat < s := by
  have := Int.emod_lt_of_pos a (by exact_mod_cast hs : (0 : ℤ) < s)
  have := Int.emod_nonneg a (by exact_mod_cast hs.ne' : (s : ℤ) ≠ 0)
  omega

lemma sumN (hs : 0 < s) (hst : s ≤ t) (hx1 : 1 ≤ xθ s t α) (x : ℕ → ℝ) (X : ℝ)
    (hX : ∀ k < s, |x k| ≤ X) (i : ℕ) :
    Summable (fun h : ℤ => eEx s t α i h * x ((((i : ℤ) + h) % s).toNat)) ∧
    |∑' h : ℤ, eEx s t α i h * x ((((i : ℤ) + h) % s).toNat)| ≤ (2 : ℝ) ^ (-xθ s t α) * X := by
  have hsum := eEx_sum (s := s) (t := t) (α := α) hs hst hx1 (i : ℤ)
  have hX0 : 0 ≤ X := (abs_nonneg _).trans (hX 0 hs)
  have hb : ∀ h : ℤ, ‖eEx s t α i h * x ((((i : ℤ) + h) % s).toNat)‖ ≤ eEx s t α i h * X :=
    fun h => by
      rw [Real.norm_eq_abs, abs_mul, abs_of_nonneg (eEx_nonneg s t α _ _)]
      exact mul_le_mul_of_nonneg_left (hX _ (toNat_emod_lt hs _)) (eEx_nonneg s t α _ _)
  have hS2 : Summable (fun h : ℤ => eEx s t α i h * X) := hsum.1.mul_right X
  have hS : Summable (fun h : ℤ => eEx s t α i h * x ((((i : ℤ) + h) % s).toNat)) :=
    hS2.of_norm_bounded hb
  refine ⟨hS, ?_⟩
  have hn : Summable (fun h : ℤ => ‖eEx s t α i h * x ((((i : ℤ) + h) % s).toNat)‖) :=
    hS2.of_nonneg_of_le (fun _ => norm_nonneg _) hb
  calc |∑' h : ℤ, eEx s t α i h * x ((((i : ℤ) + h) % s).toNat)|
      = ‖∑' h : ℤ, eEx s t α i h * x ((((i : ℤ) + h) % s).toNat)‖ := (Real.norm_eq_abs _).symm
    _ ≤ ∑' h : ℤ, ‖eEx s t α i h * x ((((i : ℤ) + h) % s).toNat)‖ := norm_tsum_le_tsum_norm hn
    _ ≤ ∑' h : ℤ, eEx s t α i h * X := hn.tsum_le_tsum hb hS2
    _ = (∑' h : ℤ, eEx s t α i h) * X := tsum_mul_right
    _ ≤ (2 : ℝ) ^ (-xθ s t α) * X := mul_le_mul_of_nonneg_right hsum.2 hX0

end glob

end J6Seg

end

section
-- Source module: Solutions.J6SegPart
/-!
# Segmented inverse (S2): the block-partition constructor

`mkBlocks s t lam w`: cut `[0, s)` at `0`, at every multiple of `lam`, and after every wrap step
(`bOf (l+1) ≠ bOf l + (t - s)`); the block of `i` starts at the last cut `≤ i` (`Nat.findGreatest`)
and ends before the next cut (`Nat.find`). All computable.

* `mkBlocks_ok`: `BlocksOK s t (mkBlocks s t lam w)` for `1 ≤ lam`, `w ≤ s`;
* `mkBlocks_len`: every block has length `≤ lam`;
* `Jlist_length_le`: `|J| ≤ s`.
-/

open Finset Real

namespace J6Seg
open IntMul.HvdH IntMul.HvdH.JPf

section part
variable (s t lam : ℕ)

/-- No wrap between `l` and `l + 1`. -/
def stepOK (l : ℕ) : Prop := bOf s t ((l : ℤ) + 1) = bOf s t (l : ℤ) + ((t : ℤ) - s)

instance (l : ℕ) : Decidable (stepOK s t l) := by unfold stepOK; infer_instance

/-- Cut points. -/
def cutP (j : ℕ) : Prop := j = 0 ∨ j % lam = 0 ∨ ¬ stepOK s t (j - 1)

instance (j : ℕ) : Decidable (cutP s t lam j) := by unfold cutP; infer_instance

/-- Start of the block of `i`. -/
def stF (i : ℕ) : ℕ := Nat.findGreatest (cutP s t lam) i

lemma enF_ex (c : ℕ) : ∃ k, cutP s t lam (c + 1 + k) ∨ s ≤ c + 1 + k := ⟨s, Or.inr (by omega)⟩

/-- `len - 1` of the block starting at `c`. -/
def enF (c : ℕ) : ℕ := Nat.find (enF_ex s t lam c)

/-- The block structure. -/
def mkBlocks (w : ℕ) : Blocks :=
  ⟨s, stF s t lam, fun i => 1 + enF s t lam (stF s t lam i), w⟩

variable {s t lam}

end part

end J6Seg

end

section
-- Source module: Solutions.J6SegTab
/-!
# Segmented inverse (S2): construction of the constant tables

`mkTables P R s t α` (computable) fills the tables of `J6SegAlg` with
* `di u = expNeg P s (α²(t-s)u²)` (`≈ 2^P e^{-πα²θu²}`),
* `dd u = expPos (P+R) s (α²(t-s)u²) R` (`≈ 2^P e^{πα²θu²}`, needs `e^{πα²θu²} ≤ 2^R`),
* `gam`: `(q;q)_u` by the recurrence `Γ_{u+1} = Γ_u (1 - q^{u+1})` (`q^k` by `expNeg`),
* `gi u = 2^P / Γ_u` (`fdiv`), `ap/am c h = (-1)^h e^{-πα²(σ ± 2β_c)h} / Γ_h`,
* `cpl i h = expNeg P s² (α² K_{i,h})` (the entries of `E`, as in `Eapprox`).

`mkTables_ok`: `TabOK` with `E = 3λ + 3` and `Dm = D_{λ-1}`, for blocks of length `≤ λ`, as soon as
`200 (3λ+1) ≤ 2^P` and `D_{λ-1} ≤ 2^R`.
-/

open Finset Real

namespace J6Seg
open IntMul.HvdH IntMul.HvdH.JPf J6QLU

/-- `max(0, ≈ 2^P q^k)`, `q = e^{-2πα² t/s}`. -/
def qN (P s t α k : ℕ) : ℤ := max (expNeg P s (2 * α ^ 2 * t * k)) 0

/-- `≈ 2^P (q;q)_u`. -/
def gamI (P s t α : ℕ) : ℕ → ℤ
  | 0 => 2 ^ P
  | u + 1 => fmul P (gamI P s t α u) (2 ^ P - qN P s t α (u + 1))

/-- `≈ 2^P / (q;q)_u`. -/
def giI (P s t α u : ℕ) : ℤ := fdiv P (2 ^ P) (gamI P s t α u)

/-- The constant tables. -/
def mkTables (P R s t α : ℕ) : Tables where
  di u := expNeg P s (α ^ 2 * (t - s) * u ^ 2)
  dd u := expPos (P + R) s (α ^ 2 * (t - s) * u ^ 2) R
  gam := gamI P s t α
  gi := giI P s t α
  ap c h := (-1) ^ h * fmul P (expNeg P s (α ^ 2 * (t + 2 * bOf s t c).toNat * h)) (giI P s t α h)
  am c h := (-1) ^ h * fmul P (expNeg P s (α ^ 2 * (t - 2 * bOf s t c).toNat * h)) (giI P s t α h)
  cpl i h := expNeg P (s ^ 2) (α ^ 2 * (eKz s t i h).toNat)

end J6Seg

end

section
-- Source module: Solutions.J6SegFull
/-!
# Segmented inverse (S2): the fully instantiated algorithm

`segInvFull P R lam w s t α U = segInvC P (mkTables P R s t α) (mkBlocks s t lam w) U`: everything
(blocks, tables, elimination) is computed from `(P, R, λ, w, s, t, α)`.

`segInvFull_err`: for `s ≤ t`, `α ≥ 2`, `α²θ ≥ 3`, `1 ≤ λ`, `w ≤ s`, `D_{λ-1} ≤ 2^R`, and
`1000 G ≤ 2^P` with `G = D_{λ-1} ((λ+1)(3λ+3)+1)(Bu+1)(λ+1)(2w+1)(s+1)²`, the output is within
`60 G/2^P + 12 e^{-πα² w(w+1)} (Bu+1) + (8/7) ε_in` of `2^P w` in both parts (twice that in norm),
where `w + Eexact w = v`.
-/

open Finset Real

namespace J6Seg
open IntMul.HvdH IntMul.HvdH.JPf J6QLU

/-- The fully instantiated segmented inverse. -/
def segInvFull (P R lam w s t α : ℕ) (U : ℕ → ℤ × ℤ) : ℕ → ℤ × ℤ :=
  segInvC P (mkTables P R s t α) (mkBlocks s t lam w) U

/-- The precision budget with `|J| ≤ s`, `E = 3λ+3`, `Dm = D_{λ-1}`. -/
noncomputable def fullG (lam w s t α : ℕ) (Bu : ℝ) : ℝ :=
  ddR (α : ℝ) ((t : ℝ) / s) (lam - 1) * (((lam : ℝ) + 1) * (3 * lam + 3) + 1) * (Bu + 1) *
    ((lam : ℝ) + 1) * (2 * w + 1) * ((s : ℝ) + 1) ^ 2

end J6Seg

end

section
-- Source module: Solutions.J6SegBand
/-!
# Segmented inverse (S2): sparse (banded + border) elimination

`geSolveS P S L A b` is the elimination of `geSolve` restricted to a pattern `S`: the Schur
update of `(i, k)` is performed only when `S i j ∧ S j k`, the right-hand side update only when
`S i j`, and the back substitution sums only over `S j k`.

`geSolveS_eq`: if `A` vanishes off `S` on `L × L`, `S` is closed under elimination in the order of
`L` (`PatOK`: `S i j → S j k → S i k` for `i, k` after the pivot `j`), and the hypotheses of
`geSolve_resid` hold (row deviation `η`, `η + |L|²/2^P ≤ 1/2`, so every pivot is positive), then
`geSolveS = geSolve`. Hence the error theorems of `geSolve` hold verbatim for the sparse version.

`bandS m c a b`: `|a - b| ≤ m` or `a ≥ c` or `b ≥ c` (band plus a border at the end, the
cyclic block-tridiagonal shape); `PatOK` holds for it along any increasing list (`bandS_patOK`).
-/

open Finset

namespace J6Seg

section band
variable (P : ℕ)

/-- Sparse elimination on a pattern `S`. -/
def geSolveS (S : ℕ → ℕ → Prop) [DecidableRel S] : List ℕ → (ℕ → ℕ → ℤ) → (ℕ → ℤ) → ℕ → ℤ
  | [], _, _ => fun _ => 0
  | j :: L, A, b =>
    let f : ℕ → ℤ := fun i => fdiv P (A i j) (A j j)
    let z' := geSolveS S L
      (fun i k => if S i j ∧ S j k then A i k - fmul P (f i) (A j k) else A i k)
      (fun i => if S i j then b i - fmul P (f i) (b j) else b i)
    fun k => if k = j then
        fdiv P (b j - ∑ k' ∈ L.toFinset.filter (fun k' => S j k'), fmul P (A j k') (z' k')) (A j j)
      else z' k

/-- Closure of a pattern under elimination in the order of `L`. -/
def PatOK (S : ℕ → ℕ → Prop) : List ℕ → Prop
  | [] => True
  | j :: L => (∀ i ∈ L, ∀ k ∈ L, S i j → S j k → S i k) ∧ PatOK S L

/-- Band of half-width `m` plus a border `≥ c`. -/
def bandS (m c : ℕ) (a b : ℕ) : Prop := (a ≤ b + m ∧ b ≤ a + m) ∨ c ≤ a ∨ c ≤ b

instance (m c : ℕ) : DecidableRel (bandS m c) := fun a b => by unfold bandS; infer_instance

end band

end J6Seg

end

section
-- Source module: Solutions.J6SegSparse
/-!
# Segmented inverse (S2): the sparse algorithm

`segInvS P lam T B U` is `segInv` with the capacitance solved by `geSolveS` on the pattern
`bandS (λ + w) (s - λ - w)`: a band of half-width `λ + w` (blocks of length `≤ λ`, coupling width
`w`) plus the border of positions `≥ s - λ - w` (the cyclic wrap). Per row of the capacitance
at most `2(λ+w)+1` band entries plus the border are touched.

* `capI_band`: the capacitance vanishes off the pattern;
* `segInvS_eq`: under the hypotheses of `segInv_err`, `segInvS = segInv`;
* `segInvFullS_err`: the fully instantiated sparse algorithm satisfies the bound of
  `segInvFull_err`.
-/

open Finset Real

namespace J6Seg
open IntMul.HvdH IntMul.HvdH.JPf J6QLU

/-- The sparse segmented inverse. -/
def segInvS (P lam : ℕ) (T : Tables) (B : Blocks) (U : ℕ → ℤ) : ℕ → ℤ :=
  let Y := applyMinv P T B U
  let z := geSolveS P (bandS (lam + B.w) (B.s - lam - B.w)) (Jlist B) (capI P T B) Y
  let Y₂ := applyMinv P T B (coupleI P T B z)
  fun i => Y i - Y₂ i

/-- Complex version. -/
def segInvCS (P lam : ℕ) (T : Tables) (B : Blocks) (U : ℕ → ℤ × ℤ) : ℕ → ℤ × ℤ :=
  fun i => (segInvS P lam T B (fun j => (U j).1) i, segInvS P lam T B (fun j => (U j).2) i)

/-- Fully instantiated sparse version. -/
def segInvFullS (P R lam w s t α : ℕ) (U : ℕ → ℤ × ℤ) : ℕ → ℤ × ℤ :=
  segInvCS P lam (mkTables P R s t α) (mkBlocks s t lam w) U

end J6Seg

end

section
-- Source module: Solutions.J6SegLn1
/-!
# Segmented inverse (S2), machine side 1: the block solve as list functions

The block solve `blkI` (seven fixed-point stages) on a list `x` of length `n`, with the tables as
lists (`di, gi, gam, dd` of length `≥ n`, the block's `am, ap` of length `≥ n`):

* `dgL c x = zipWith (fmul P) c x` (diagonal stage);
* `convL c x` = the exact lower-triangular Toeplitz sums `∑_{k ≤ u} c_{u-k} x_k` (`u < |x|`),
  `lwL c x = (convL c x).map (rs P)` (one rounding per output);
* `upL c x = (lwL c x.reverse).reverse` (the upper stage is the lower stage on the reversed list);
* `blkL = dgL dd ∘ dgL gi ∘ upL ap ∘ dgL gam ∘ lwL am ∘ dgL gi ∘ dgL di`.

**`blkL_getD`**: `blkL` computes `blkI` entrywise.
-/

open Finset

namespace J6Seg

section lst
variable (P : ℕ)

/-- Exact lower-triangular Toeplitz sums. -/
def convL (c x : List ℤ) : List ℤ :=
  (List.range x.length).map fun u => ∑ k ∈ range (u + 1), c.getD (u - k) 0 * x.getD k 0

/-- Diagonal stage on lists. -/
def dgL (c x : List ℤ) : List ℤ := List.zipWith (fmul P) c x

/-- Lower Toeplitz stage on lists. -/
def lwL (c x : List ℤ) : List ℤ := (convL c x).map (rs P)

/-- Upper Toeplitz stage on lists. -/
def upL (c x : List ℤ) : List ℤ := (lwL P c x.reverse).reverse

/-- The block solve on lists. -/
def blkL (di gi gam dd am ap x : List ℤ) : List ℤ :=
  dgL P dd (dgL P gi (upL P ap (dgL P gam (lwL P am (dgL P gi (dgL P di x))))))

/-- A list represents a function on `[0, n)`. -/
def RepL (n : ℕ) (L : List ℤ) (g : ℕ → ℤ) : Prop := L.length = n ∧ ∀ u < n, L.getD u 0 = g u

end lst

end J6Seg

end

section
-- Source module: Solutions.J6SegLn2
/-!
# Segmented inverse (S2), machine side 2: `M⁻¹` chunk by chunk

With `λ`-aligned chunks `[cλ, cλ + λ)` every block of `mkBlocks s t λ w` lies inside one chunk
(cuts at the multiples of `λ`). The blocks of chunk `c` (`chunkBl`) start at `cλ` and follow each
other (`b ↦ b + 1 + enF b`) until `min(s, cλ + λ)`.

* `chunkL di gi gam dd descs x`: the block solves of consecutive blocks
  `descs = [(n_k, am_k, ap_k)]` on consecutive segments of `x`, zeros on the rest;
* `descsOf T c`: the block descriptors of chunk `c` (lengths and the block tables `am, ap`);
* **`chunkL_rep`**: on chunk `c`, `chunkL … (descsOf c) x` is `M⁻¹` (`applyMinv`) at the positions
  `< s` of the chunk and `0` at the others.
-/

open Finset

namespace J6Seg

section chunk
variable (P : ℕ) (s t lam : ℕ)

/-- Blocks `(start, length)` from `st` up to the end `e` (fuel `f`). -/
def blFrom (e : ℕ) : ℕ → ℕ → List (ℕ × ℕ)
  | 0, _ => []
  | f + 1, st => if st < e then (st, 1 + enF s t lam st) :: blFrom e f (st + 1 + enF s t lam st)
      else []

/-- The blocks of chunk `c`. -/
def chunkBl (c : ℕ) : List (ℕ × ℕ) := blFrom s t lam (min s (c * lam + lam)) lam (c * lam)

/-- The block descriptors of chunk `c`: length and the tables `am, ap` of the block. -/
def descsOf (T : Tables) (c : ℕ) : List (ℕ × List ℤ × List ℤ) :=
  (chunkBl s t lam c).map fun b => (b.2, (List.range b.2).map (T.am b.1),
    (List.range b.2).map (T.ap b.1))

/-- The chunk kernel on lists. -/
def chunkL (di gi gam dd : List ℤ) : List (ℕ × List ℤ × List ℤ) → List ℤ → List ℤ
  | [], x => x.map fun _ => 0
  | (n, am, ap) :: ds, x => blkL P di gi gam dd am ap (x.take n) ++ chunkL di gi gam dd ds (x.drop n)

variable {s t lam}

end chunk

section chunkRep
variable (P R s t lam w α : ℕ)

/-- The tables and blocks of the line. -/
abbrev TT : Tables := mkTables P R s t α
abbrev BB : Blocks := mkBlocks s t lam w

/-- The global stage tables as lists of length `λ`. -/
def diL : List ℤ := (List.range lam).map (TT P R s t α).di
def giL : List ℤ := (List.range lam).map (TT P R s t α).gi
def gamL : List ℤ := (List.range lam).map (TT P R s t α).gam
def ddL : List ℤ := (List.range lam).map (TT P R s t α).dd

/-- The chunk kernel with the line's tables. -/
def chunkK (c : ℕ) (x : List ℤ) : List ℤ :=
  chunkL P (diL P R s t lam α) (giL P R s t lam α) (gamL P R s t lam α) (ddL P R s t lam α)
    (descsOf s t lam (TT P R s t α) c) x

/-- The descriptors of the blocks from `st`. -/
def descsFrom (c f st : ℕ) : List (ℕ × List ℤ × List ℤ) :=
  (blFrom s t lam (min s (c * lam + lam)) f st).map fun b =>
    (b.2, (List.range b.2).map ((TT P R s t α).am b.1), (List.range b.2).map ((TT P R s t α).ap b.1))

variable {P R s t lam w α}

end chunkRep

end J6Seg

end

section
-- Source module: Solutions.J6SegLn3
/-!
# Segmented inverse (S2), machine side 3: the sparse elimination in closed form

For the elimination `geSolveS P S L A b` along `L = [l₀, l₁, …]` (no repetitions), with the
matrices `Aseq L A r` and right-hand sides `bseq L A b r` after `r` pivots:

* **`bseq_closed`**: `bseq r i = b i − ∑_{q < r, S i l_q} fmul (F_q i) (bt q)` where
  `F_q i = fdiv (A_q i l_q) (A_q l_q l_q)` and `bt q = bseq q l_q` (forward substitution: every
  pivot step subtracts one rounded term);
* **`geSolveS_closed`**: `z l_r = fdiv (bt r − ∑_{k ∈ tail, S l_r k} fmul (A_r l_r k) (z k)) (A_r l_r l_r)`
  (back substitution with the final values `z`);
* `geSolveS_off`: `z k = 0` off `L`.

Everything data-independent (`A_q`, hence all multipliers and pivots) can be precomputed; the
data enter only through `bt` and `z`.
-/

open Finset

namespace J6Seg

section closed
variable (P : ℕ) (S : ℕ → ℕ → Prop) [DecidableRel S]

/-- One Schur step of the pattern-restricted elimination. -/
def stepA (j : ℕ) (A : ℕ → ℕ → ℤ) : ℕ → ℕ → ℤ :=
  fun i k => if S i j ∧ S j k then A i k - fmul P (fdiv P (A i j) (A j j)) (A j k) else A i k

/-- One right-hand-side step. -/
def stepB (j : ℕ) (A : ℕ → ℕ → ℤ) (b : ℕ → ℤ) : ℕ → ℤ :=
  fun i => if S i j then b i - fmul P (fdiv P (A i j) (A j j)) (b j) else b i

/-- The matrix after `r` pivots of `L`. -/
def Aseq : List ℕ → (ℕ → ℕ → ℤ) → ℕ → ℕ → ℕ → ℤ
  | [], A, _ => A
  | _ :: _, A, 0 => A
  | j :: L, A, r + 1 => Aseq L (stepA P S j A) r

/-- The right-hand side after `r` pivots of `L`. -/
def bseq : List ℕ → (ℕ → ℕ → ℤ) → (ℕ → ℤ) → ℕ → ℕ → ℤ
  | [], _, b, _ => b
  | _ :: _, _, b, 0 => b
  | j :: L, A, b, r + 1 => bseq L (stepA P S j A) (stepB P S j A b) r

end closed

end J6Seg

end

section
-- Source module: Solutions.J6SegLn4
/-!
# Segmented inverse (S2), machine side 4: the capacitance solve as list kernels

The per-line work on the near-boundary positions `L = Jlist B = [j₀ < j₁ < …]` is done by two
list kernels driven by precomputed (data-independent) coefficient streams:

* `fwdK fs Y []`: for `r = 0, 1, …`: `bt_r = Y_r − ∑_i fmul (fs_r)_i bt_{r-1-i}` (the accumulator
  holds the previous values, most recent first); the result is `[bt_{m-1}, …, bt_0]`;
* `bwdK bs bt [] []`: for `r = m-1, …, 0` with the record `(cb, cz, pv, fl)`:
  `z_r = fdiv (bt_r − ∑ fmul cb_i z_{r+1+i} − ∑ fmul cz_i ZB_i) pv`, where `ZB` holds the
  `z` of the border indices computed so far (`fl = 1` marks a border index); the result is
  `[z₀, …, z_{m-1}]`.

**`fwdK_eq`**, **`bwdK_eq`**: with the coefficient streams `fsOf`, `bsOf` of the line (built from
`Aseq` of the pattern-restricted elimination of `capI`), `fwdK` computes the forward values
`bseq r j_r` and `bwdK` the solution `geSolveS … j_r`, given the support facts of the pattern
(`FwdSupp`, `BwdSupp`).
-/

open Finset

namespace J6Seg

section kern
variable (P : ℕ)

/-- `∑_i fmul cs_i acc_i`. -/
def dotF (cs acc : List ℤ) : ℤ := ∑ i ∈ range cs.length, fmul P (cs.getD i 0) (acc.getD i 0)

/-- The forward kernel. -/
def fwdK : List (List ℤ) → List ℤ → List ℤ → List ℤ
  | cs :: css, y :: ys, acc => fwdK css ys ((y - dotF P cs acc) :: acc)
  | _, _, acc => acc

/-- The backward kernel (records `(cb, cz, pv, fl)`). -/
def bwdK : List (List ℤ × List ℤ × ℤ × ℤ) → List ℤ → List ℤ → List ℤ → List ℤ
  | (cb, cz, pv, fl) :: rs, bt :: bts, accZ, accB =>
    bwdK rs bts (fdiv P (bt - dotF P cb accZ - dotF P cz accB) pv :: accZ)
      (if fl = 1 then fdiv P (bt - dotF P cb accZ - dotF P cz accB) pv :: accB else accB)
  | _, _, accZ, _ => accZ

end kern

section fwd
variable (P : ℕ) (S : ℕ → ℕ → Prop) [DecidableRel S] (L : List ℕ) (A : ℕ → ℕ → ℤ) (b : ℕ → ℤ)

/-- The forward values `bt_r = bseq r l_r`. -/
def btV (r : ℕ) : ℤ := bseq P S L A b r (L.getD r 0)

/-- The forward multiplier of row `r` on pivot `q`. -/
def fco (r q : ℕ) : ℤ :=
  if S (L.getD r 0) (L.getD q 0) then
    fdiv P (Aseq P S L A q (L.getD r 0) (L.getD q 0)) (Aseq P S L A q (L.getD q 0) (L.getD q 0))
  else 0

/-- The forward coefficient list of row `r` with window `c` (most recent pivot first). -/
def fsW (r c : ℕ) : List ℤ := (List.range c).map fun i => fco P S L A r (r - 1 - i)

variable {P S L A b}

/-- Support of a window: no pattern entry before `r - c`. -/
def FwdSupp (r c : ℕ) : Prop := c ≤ r ∧ ∀ q < r - c, ¬ S (L.getD r 0) (L.getD q 0)

/-- The accumulator after `k` forward rows. -/
def accF (k : ℕ) : List ℤ := ((List.range k).map (btV P S L A b)).reverse

end fwd

section bwd
variable (P : ℕ) (S : ℕ → ℕ → Prop) [DecidableRel S] (L : List ℕ) (A : ℕ → ℕ → ℤ) (b : ℕ → ℤ)

/-- The solution at position `r` of the line. -/
def zW (r : ℕ) : ℤ := geSolveS P S L A b (L.getD r 0)

/-- `[z_k, …, z_{m-1}]`. -/
def accZ (k : ℕ) : List ℤ := (List.range' k (L.length - k)).map (zW P S L A b)

/-- The border values among `z_k, …, z_{m-1}` (in increasing order). -/
def accB (fl : ℕ → ℤ) (k : ℕ) : List ℤ :=
  ((List.range' k (L.length - k)).filter (fun q => fl q = 1)).map (zW P S L A b)

/-- The back-substitution sum of row `r`. -/
def bsum (r : ℕ) : ℤ :=
  ∑ k' ∈ (L.drop (r + 1)).toFinset.filter (fun k' => S (L.getD r 0) k'),
    fmul P (Aseq P S L A r (L.getD r 0) k') (geSolveS P S L A b k')

variable {P S L A b}

end bwd

end J6Seg

end

section
-- Source module: Solutions.J6SegLn5
/-!
# Segmented inverse (S2), machine side 5: the coefficient streams of the backward kernel

For a sorted line `L` and a border predicate `bd` (border positions: the cyclic wrap), the
back-substitution sum of row `r` splits into

* a band part over the next `W` rows (`cbOf W r`, zero on border rows and off the pattern), and
* a border part over all later border rows (`czOf r`, against `accB` with flags `flOf`).

**`hsplit_of`**: the split is exact (the hypothesis of `bwdK_eq`), given the window property
`hwin`: a non-border pattern entry `(L_r, L_q)` with `r < q` has `q ≤ r + W`.
**`win_of_band`**: for a strictly increasing `L`, a band of half-width `W` in positions gives
the window `W` in indices.
-/

open Finset

namespace J6Seg

section split
variable (P : ℕ) (S : ℕ → ℕ → Prop) [DecidableRel S] (L : List ℕ) (A : ℕ → ℕ → ℤ)
  (bd : ℕ → Prop) [DecidablePred bd]

/-- Border flags. -/
def flOf (q : ℕ) : ℤ := if bd (L.getD q 0) then 1 else 0

/-- Band coefficients of row `r` (window `W`). -/
def cbOf (W r : ℕ) : List ℤ := (List.range W).map fun i =>
  if r + 1 + i < L.length ∧ ¬ bd (L.getD (r + 1 + i) 0) ∧ S (L.getD r 0) (L.getD (r + 1 + i) 0)
  then Aseq P S L A r (L.getD r 0) (L.getD (r + 1 + i) 0) else 0

/-- The later border indices of row `r`. -/
def bIdx (r : ℕ) : List ℕ :=
  (List.range' (r + 1) (L.length - (r + 1))).filter (fun q => flOf L bd q = 1)

/-- Border coefficients of row `r`. -/
def czOf (r : ℕ) : List ℤ := (bIdx L bd r).map fun q =>
  if S (L.getD r 0) (L.getD q 0) then Aseq P S L A r (L.getD r 0) (L.getD q 0) else 0

variable {P S L A bd}

end split

end J6Seg

end

section
-- Source module: Solutions.J6SegLn6
/-!
# Segmented inverse (S2), machine side 6: the coupling as a cyclic window

On a sorted line `L = [j₀ < … < j_{m-1}]` of positions `< s` (the near-boundary positions), with
`z` vanishing off `L`, the coupling `W = U z` at `j_r` only reads `z` at the line indices
`q` within cyclic index distance `K ≥ w` of `r` (positions within cyclic distance `w`, and the
line is strictly increasing). With the cyclically extended list
`zext = (last K of zl) ++ zl ++ (first K of zl)`, `zl = [z_{j_0}, …, z_{j_{m-1}}]`:

* **`couple_window`**: `W_{j_r} = ∑_{i ≤ 2K} fmul (cc_r)_i zext_{r+i}` (`2K+1 ≤ m`, `2w < s`);
* **`couple_off`**: `W_i = 0` at positions that are not near a block boundary.
-/

open Finset

namespace J6Seg

section cpl
variable (P : ℕ) (T : Tables) (B : Blocks) (L : List ℕ) (K : ℕ)

/-- The values on the line. -/
def zlOf (zf : ℕ → ℤ) : List ℤ := (List.range L.length).map fun q => zf (L.getD q 0)

/-- The cyclic extension. -/
def zext (zl : List ℤ) : List ℤ := zl.drop (zl.length - K) ++ (zl ++ zl.take K)

/-- The line index read by window slot `i` of row `r`. -/
def qOf (r i : ℕ) : ℕ := (r + i + L.length - K) % L.length

/-- The window coefficients of row `r`. -/
def ccOf (r i : ℕ) : ℤ := ∑ h ∈ Hset B (L.getD r 0),
  if col B.s (L.getD r 0) h = L.getD (qOf L K r i) 0 then T.cpl (L.getD r 0) h else 0

variable {P T B L K}

end cpl

end J6Seg

end

section
-- Source module: Solutions.J6SegLn7
/-!
# Segmented inverse (S2), machine side 7: the capacitance solve of one line as list kernels

For the pattern `S = bandS (λ+w) c` (`c = s - λ - w`, border `bd = (c ≤ ·)`) on the sorted
line `L = Jlist B`:

* `WbOf r` (backward band window): the number of later line positions within `λ + w` of `L_r`;
* `cfOf r` (forward window): all previous rows on border rows, else the number of previous line
  positions within `λ + w`;
* **`jsolve_eq`**: `bwdK (records) (fwdK (forward streams) [Y_{L_r}]_r [])` is the list of the
  solution `geSolveS … Y` at `L_0, …, L_{m-1}`, for every right-hand side `Y`.

All windows and streams depend only on `s, t, λ, w` and the tables (not on the data).
-/

open Finset

namespace J6Seg

section jsolve
variable (P : ℕ) (L : List ℕ) (mw c : ℕ) (A : ℕ → ℕ → ℤ)

/-- Backward band window of row `r`. -/
def WbOf (r : ℕ) : ℕ := ((Ioo r L.length).filter (fun q => L.getD q 0 ≤ L.getD r 0 + mw)).card

/-- Forward window of row `r`. -/
def cfOf (r : ℕ) : ℕ :=
  if c ≤ L.getD r 0 then r else ((range r).filter (fun q => L.getD r 0 ≤ L.getD q 0 + mw)).card

/-- The forward coefficient streams. -/
def fsOf : List (List ℤ) :=
  (List.range L.length).map fun r => fsW P (bandS mw c) L A r (cfOf L mw c r)

/-- The backward records (rows in decreasing order). -/
def bsOf : List (List ℤ × List ℤ × ℤ × ℤ) :=
  (List.range L.length).reverse.map fun r =>
    (cbOf P (bandS mw c) L A (c ≤ ·) (WbOf L mw r) r, czOf P (bandS mw c) L A (c ≤ ·) r,
      Aseq P (bandS mw c) L A r (L.getD r 0) (L.getD r 0), flOf L (c ≤ ·) r)

variable {P L mw c A}

end jsolve

end J6Seg

end

section
-- Source module: Solutions.J6SegLn8
/-!
# Segmented inverse (S2), machine side 8: one real line of `segInvS` by list kernels

The real part of the segmented inverse on one line, as the machine computes it, on a padded
line `x` of `N = n_c λ` positions (`s ≤ N`):

1. `Y = applyL x`: the chunk kernel on each `λ`-chunk (`M⁻¹`);
2. the values of `Y` on the line `L = Jlist B`, the two capacitance kernels (`jsolve_eq`);
3. the coupling window kernel on the cyclic extension (`couple_window`), scattered back to the
   positions (zero off `L`);
4. `Y₂ = applyL` of the scattered coupling; the output is `Y - Y₂`.

**`segLineC_eq`**: entry `u < N` of the result is `segInvS … u` for `u < s` and `0` beyond.
-/

open Finset

namespace J6Seg

section chunks

/-- A kernel applied to consecutive `λ`-chunks. -/
def chunksL (f : ℕ → List ℤ → List ℤ) (lam nc : ℕ) (x : List ℤ) : List ℤ :=
  (List.range nc).flatMap fun c => f c ((x.drop (c * lam)).take lam)

end chunks

section line
variable (P R s t lam w α : ℕ)

/-- `M⁻¹` on a padded line by chunks. -/
def applyL (nc : ℕ) (x : List ℤ) : List ℤ := chunksL (chunkK P R s t lam α) lam nc x

variable {P R s t lam w α}

end line

section segline
variable (P R s t lam w α K : ℕ)

/-- The line of near-boundary positions. -/
abbrev LJ : List ℕ := Jlist (BB s t lam w)

/-- The coupling window coefficients of row `r`. -/
def ccL (r : ℕ) : List ℤ := (List.range (2 * K + 1)).map (ccOf (TT P R s t α) (BB s t lam w)
  (LJ s t lam w) K r)

/-- **One real line of the segmented inverse, by list kernels.** -/
def segLineC (nc : ℕ) (x : List ℤ) : List ℤ :=
  let L := LJ s t lam w
  let A := capI P (TT P R s t α) (BB s t lam w)
  let Yl := applyL P R s t lam α nc x
  let zl := bwdK P (bsOf P L (lam + w) (s - lam - w) A)
    (fwdK P (fsOf P L (lam + w) (s - lam - w) A)
      ((List.range L.length).map fun r => Yl.getD (L.getD r 0) 0) []) [] []
  let Wl := (List.range L.length).map fun r => dotF P (ccL P R s t lam w α K r) ((zext K zl).drop r)
  let Wp := (List.range (nc * lam)).map fun u => if u ∈ L then Wl.getD (L.idxOf u) 0 else 0
  let Y2 := applyL P R s t lam α nc Wp
  (List.range (nc * lam)).map fun u => Yl.getD u 0 - Y2.getD u 0

/-- The first solve on a line. -/
def YsolL (x : List ℤ) : ℕ → ℤ := applyMinv P (TT P R s t α) (BB s t lam w) (fun v => x.getD v 0)

/-- The capacitance solution on a line. -/
def zsolL (x : List ℤ) : ℕ → ℤ :=
  geSolveS P (bandS (lam + w) (s - lam - w)) (LJ s t lam w) (capI P (TT P R s t α) (BB s t lam w))
    (YsolL P R s t lam w α x)

variable {P R s t lam w α K}

end segline

end J6Seg

end

section
-- Source module: Solutions.J6SegB0
/-!
# Segmented inverse (S2): Proposition 4.7 and Lemma 19 with the segmented `J̃`

The base `Bapprox = P̃_s⁻¹ D̃' J̃' C̃ P̃_t` (HvdH Proposition 4.7) solves the resampling system
`(I + E) w = y/2` by the Neumann iteration `Japprox`. Here the solver is the segmented inverse
`segInvFullS` (sparse banded elimination, tables and blocks computed from the parameters):

* `segJ p P R λ w s t α Y`: scale `Y` (precision `p`) to precision `P` (exact shift by
  `P - p - 1`, which also halves), run `segInvFullS`, shift back by `P - p` (floor), and read
  `0` outside `[0, s)`.
* `Bseg = P̃_s⁻¹ D̃' segJ C̃ P̃_t` and the permutation-free `B0seg = D̃' segJ C̃` (`B0segP`: read
  mod `s`).
* `prop_4_7_seg`: the maps `A`, `B` of Proposition 4.7 (built exactly as in the base proof), with
  `‖toC (Bseg V)_j - 2^p (B v)_j‖ ≤ 2 segE + 3`, where
  `segE = 2^p (60 G/2^P + 18 e^{-πα² w(w+1)}) + 1`, `G = fullG λ w s t α (1/2)`.
  Hypotheses: `α²θ ≥ 3` (`3 s ≤ α² (t - s)`), `α² ≤ p`, `1 ≤ λ`, `w ≤ s`, `D_{λ-1} ≤ 2^R`,
  `p + 1 ≤ P`, `1000 G ≤ 2^P`, `200 segE ≤ 2^p` (no `jIter` hypothesis).
* `lemma19_seg`: the same in the `B₀ = B0of s t B` form with the `permS`/`permT` identity of
  `AxisMaps6.hid` / `KFM.AxisMaps.hid`.
-/

open Complex Real

namespace J6Seg
open IntMul.HvdH IntMul.HvdH.Res IntMul.HvdH.P47

/-- `segJ`: the segmented solver on `p`-bit data (`P`-bit internal precision). -/
def segJ (p P R lam w s t α : ℕ) (Y : ℕ → ℤ × ℤ) (l : ℕ) : ℤ × ℤ :=
  if l < s then
    ((segInvFullS P R lam w s t α
        (fun i => ((Y i).1 * 2 ^ (P - p - 1), (Y i).2 * 2 ^ (P - p - 1))) l).1 / 2 ^ (P - p),
     (segInvFullS P R lam w s t α
        (fun i => ((Y i).1 * 2 ^ (P - p - 1), (Y i).2 * 2 ^ (P - p - 1))) l).2 / 2 ^ (P - p))
  else (0, 0)

/-- `B̃ = P̃_s⁻¹ D̃' segJ C̃ P̃_t`. -/
def Bseg (p P R lam w s t α : ℕ) (V : ℕ → ℤ × ℤ) : ℕ → ℤ × ℤ :=
  psInvPerm s t (Dapprox p s t α (segJ p P R lam w s t α (cSel s t (ptPerm s t V))))

/-- The permutation-free `B̃₀ = D̃' segJ C̃`. -/
def B0seg (p P R lam w s t α : ℕ) (W : ℕ → ℤ × ℤ) : ℕ → ℤ × ℤ :=
  Dapprox p s t α (segJ p P R lam w s t α (cSel s t W))

/-- The error of `segJ` per real part, in units of `2^{-p}`. -/
noncomputable def segE (p P lam w s t α : ℕ) : ℝ :=
  2 ^ p * (60 * fullG lam w s t α (1 / 2) / 2 ^ P +
    3 * (4 * Real.exp (-(π * (α : ℝ) ^ 2 * (w * (w + 1))))) * (1 / 2 + 1)) + 1

/-! ### Lemma 19 form (`B₀ = P_s B P_t⁻¹`) -/

end J6Seg

end

section
-- Source module: Solutions.J6SegLn9
/-!
# Segmented inverse (S2), machine side 9: one complex line of `B₀`

The whole `B₀` line map (`B0seg = D̃' segJ C̃`) on a line `y` of `t` Gaussian items, as list
passes:

1. selection and scaling (`selC`): row `l < s` is the item `⌊(2tl + s)/(2s)⌋` times
   `2^{P-p-1}`, rows `s ≤ l < N` are `0`;
2. `segLineC` on the real and on the imaginary parts;
3. the output pass: floor by `2^{P-p}`, multiply by `d̃'_l`, round to zero by `2^p`
   (rows `≥ s` are `0`).

**`lineB_eq`**: entry `l < N` of `lineB y` is `B0seg … (y as a function) l`.
-/

open Finset

namespace J6Seg

section lineB
variable (p P R s t lam w α K : ℕ)

/-- The selected row of `l`. -/
def rowSel (l : ℕ) : ℕ := (2 * t * l + s) / (2 * s)

/-- Pass 1 (one part): select and scale. -/
def selC (N : ℕ) (g : ℤ × ℤ → ℤ) (y : List (ℤ × ℤ)) : List ℤ :=
  (List.range N).map fun l => if l < s then g (y.getD (rowSel s t l) (0, 0)) * 2 ^ (P - p - 1) else 0

/-- The output pass on one entry. -/
def outE (l : ℕ) (x : ℤ) : ℤ := (IntMul.HvdH.dTilde p s t α l * (x / 2 ^ (P - p))).tdiv (2 ^ p)

/-- **The `B₀` line map by list passes** (`N = n_c λ` rows). -/
def lineB (nc : ℕ) (y : List (ℤ × ℤ)) : List (ℤ × ℤ) :=
  let xr := segLineC P R s t lam w α K nc (selC p P s t (nc * lam) Prod.fst y)
  let xi := segLineC P R s t lam w α K nc (selC p P s t (nc * lam) Prod.snd y)
  (List.range (nc * lam)).map fun l =>
    if l < s then (outE p P s t α l (xr.getD l 0), outE p P s t α l (xi.getD l 0)) else (0, 0)

variable {p P R s t lam w α K}

end lineB

end J6Seg

end

section
-- Source module: Solutions.J6SegLn11
/-!
# Segmented inverse (S2), machine side 11: the `B₀` tensor pass

Axis by axis (`tensorGo`, as the frozen `KFY.passB`), every `t_i`-line `X` is replaced by the
`t_i` values `B0seg … (X as a function) l` (zero for `l ≥ s_i`): this is what `lineB` computes
(`lineB_list`, with `n_c λ = t_i`).

**`passBseg`**: for per-axis maps `Bm i` with `‖Bm i‖ ≤ 1/2`, the per-axis bound `9` of
`r6_axisMaps_fit` and `B0seg 0 = 0`, from `RepFull ys y`, `‖y‖ ≤ 1`:
`RepPad` of the result, every entry in the disc (`GB (p+2)`), and
`‖vOf − tensCLM (kerB Bm) y‖ ≤ 9 d · 2^{-p}`.
-/

set_option linter.unusedSectionVars false
set_option linter.unusedVariables false

open Complex

namespace J6Seg
open IntMul IntMul.HvdH IntMul.HvdH.Ref IntMul.HvdH.RefPf NLogN.KFM NLogN.KFY IntMul.OpenAI.RT
open NLogN.KFB (GI gflat)

/-! ## `lineB` as a list map -/

/-! ## Facts on `B0seg` -/

/-! ## Padded vs periodic slice maps -/

/-! ## The tensor pass -/

section pass
variable {d : ℕ} (p R lam w α : ℕ) (Pf : ℕ → ℕ) (s t : Fin d → ℕ) [∀ i, NeZero (s i)]
  [∀ i, NeZero (t i)]

/-- The `B₀` map of axis `i` (as a function of the axis index in `ℕ`). -/
def B0ax (i : ℕ) : (ℕ → ℤ × ℤ) → (ℕ → ℤ × ℤ) :=
  B0seg p (Pf ((List.ofFn s).getD i 1)) R lam w ((List.ofFn s).getD i 1)
    ((List.ofFn t).getD i 1) α

/-- The padded maps on `Fin d`. -/
def B0F : Fin d → (ℕ → ℤ × ℤ) → (ℕ → ℤ × ℤ) := fun i => B0ax p R lam w α Pf s t i

/-- The periodic maps on `Fin d` (read mod `s_i`). -/
def B0G : Fin d → (ℕ → ℤ × ℤ) → (ℕ → ℤ × ℤ) := fun i U l => B0ax p R lam w α Pf s t i U (l % s i)

/-- The `B₀` pass maps for `tensorGo`. -/
def mapsBs : List (ℕ × (List GI → List GI)) :=
  (List.range d).map fun i => ((List.ofFn t).getD i 1, fun X : List GI =>
    (List.range ((List.ofFn t).getD i 1)).map (B0ax p R lam w α Pf s t i fun j => X.getD j (0, 0)))

end pass

end J6Seg

end

section
-- Source module: Solutions.J6SegLn10
/-!
# Segmented inverse (S2), machine side 10: the solver maps zero to zero

`B₀` must send an all-zero line to an all-zero line (the padded entries of the tensor stay
exactly zero). The only non-obvious step is the rounded division `fdiv P 0 d = 0`, true for
`d ≥ 0`: every pivot of the elimination is positive in the near-identity regime.

* `PivPos`, **`geSolveS_zero`**: nonnegative pivots give `geSolveS … 0 = 0`;
* **`pivPos_of`**: the hypotheses of `geSolveS_eq` give nonnegative pivots (same induction);
* **`segInvS_zero`** (hypotheses of `segInvS_eq`), **`segInvFullS_zero`** (parameter
  hypotheses of `segInvFullS_err`).
-/

open Finset Real

namespace J6Seg
open IntMul.HvdH IntMul.HvdH.JPf J6QLU

section piv
variable (P : ℕ) (S : ℕ → ℕ → Prop) [DecidableRel S]

/-- Nonnegative pivots along the elimination order. -/
def PivPos : List ℕ → (ℕ → ℕ → ℤ) → Prop
  | [], _ => True
  | j :: L, A => 0 ≤ A j j ∧ PivPos L (stepA P S j A)

end piv

end J6Seg

end

section
-- Source module: Solutions.J6SegFitA
/-!
# Segmented inverse (S2): fitted parameters and the `D_{λ-1} ≤ 2^R` bound

`lemma19_seg_fit`: any `1 ≤ λ ≤ d k`, `w = 7k` (`k = ⌊√d⌋ + 1`), `R = ⌊p/7⌋`, `P = segPF p d s`,
error `≤ 9`, `‖B₀‖ ≤ 1/2`, and `B₀ 0 = 0` exactly (`B0seg_zero`). `lamP d = 2^{⌊log₂ (d k)⌋}` is the power-of-two choice used by the
machine (blocks never cross the `λ`-aligned chunks).
-/

open Complex Real

namespace J6Seg
open IntMul.HvdH IntMul.HvdH.Res IntMul.HvdH.P47
/-! ### Fitted parameters (cost row `d p^δ`): `λ = d k`, `w = 7k`, `k = ⌊√d⌋ + 1`, `R = ⌊p/7⌋` -/

/-- `k = ⌊√d⌋ + 1` (`d < k² ≤ 4d`). -/
def kF (d : ℕ) : ℕ := Nat.sqrt d + 1
/-- The sub-block length `λ = d k ≈ d^{3/2}`. -/
def lamF (d : ℕ) : ℕ := d * kF d
/-- The coupling width `w = 7k ≈ 7√d`. -/
def wF (d : ℕ) : ℕ := 7 * kF d
/-- The budget numerator `1000 G / D_{λ-1}`. -/
def NF (d s : ℕ) : ℕ :=
  1500 * ((lamF d + 1) * (3 * lamF d + 3) + 1) * (lamF d + 1) * (2 * wF d + 1) * (s + 1) ^ 2
/-- The working precision `P = p + 1 + ⌊p/7⌋ + size (NF d s)`. -/
def segPF (p d s : ℕ) : ℕ := p + 1 + p / 7 + Nat.size (NF d s)

/-- The power-of-two sub-block length `λ = 2^{⌊log₂ (d k)⌋}` (`d k / 2 < λ ≤ d k`). -/
def lamP (d : ℕ) : ℕ := 2 ^ Nat.log 2 (lamF d)

end J6Seg

end