Round-six κ development, definitions part 3

-- 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_Numeric
import Definitions.Def_IntMul_HvdH_ExpApprox
import Definitions.Def_IntMul_HvdH_Twiddle
import Definitions.Def_IntMul_HvdH_RealExp
import Definitions.Def_IntMul_HvdH_EApprox
import Definitions.Def_IntMul_HvdH_JApprox

set_option autoImplicit false
set_option relaxedAutoImplicit false
set_option maxHeartbeats 1000000

section
-- Source module: Solutions.J6QLUBasic
/-!
# q-Pochhammer symbols and Gaussian binomials (algebraic core of the segmented-inverse q-LU)

Over an arbitrary field `K`:
* `qpoch q n = (q;q)_n = ∏_{i<n} (1 - q^(i+1))`,
* `qbin q n k` = Gaussian binomial `[n k]_q`, defined by the q-Pascal recurrence
  `[n+1, k+1] = [n, k] + q^(k+1) [n, k+1]`,
* `qbin_mul_qpoch` : `[n k] (q;q)_k (q;q)_{n-k} = (q;q)_n`,
* `qbin_rothe` (q-binomial theorem / Rothe):
  `∑_{j ≤ h} (-1)^j q^(j choose 2) [h j] y^j = ∏_{i<h} (1 - y q^i)`,
* `qbin_alt_sum` : `∑_{j ≤ h} (-1)^j q^(j choose 2) [h j] = δ_{h0}`,
* `qpascal_identity` (q-Newton type identity behind the LU factor):
  `∑_{k ≤ u} [u k] (q^(v-k+1);q)_k q^((u-k)(v-k)) = 1`.
-/

open Finset

namespace J6QLU

variable {K : Type*} [Field K]

@[simp] lemma qpoch_zero (q : K) : qpoch q 0 = 1 := by simp [qpoch]

@[simp] lemma qbin_zero_right (q : K) (n : ℕ) : qbin q n 0 = 1 := by
  cases n <;> simp [qbin]

@[simp] lemma qbin_zero_succ (q : K) (k : ℕ) : qbin q 0 (k + 1) = 0 := by simp [qbin]

end J6QLU

end

section
-- Source module: Solutions.J6SegFix
/-!
# Segmented inverse (S2): fixed-point primitives and elimination without pivoting

Fixed-point numbers of precision `P` are integers `a` standing for `a / 2^P`.
* `rq a b = ⌊a/b + 1/2⌋` (rounded quotient), `rs P a = rq a (2^P)`,
  `fmul P x y = rs P (x*y)` and `fdiv P x y = rq (x 2^P) y` (errors `≤ 1/2` ulp).
* `geSolve P L A b`: Gaussian elimination **without pivoting** in the order of the index list
  `L`, written as forward elimination of the head pivot followed by back substitution.
* `schur_row` (report L9): if every row of `A - I` has absolute sum `≤ η ≤ 1/2`, the same holds for
  the Schur complement after eliminating one pivot; hence for every Schur complement, and no
  pivoting is ever needed (any order, any sparsity pattern, cyclic borders included).
* `near_id`: `‖A - I‖_∞ ≤ η < 1` and `‖A z - r‖ ≤ ρ` give `‖z‖ ≤ (‖r‖ + ρ)/(1 - η)`.
* `geSolve_resid`: the computed solution has residual `≤ (Z+1)(|L|+1)²/2^P`;
  `geSolve_bound`: and `|z| ≤ 4B + 1`.
-/

open Finset

namespace J6Seg

/-- Rounded quotient `⌊a/b + 1/2⌋` (meaningful for `b > 0`). -/
def rq (a b : ℤ) : ℤ := (2 * a + b) / (2 * b)

/-- Rounded right shift `⌊a/2^P + 1/2⌋`. -/
def rs (P : ℕ) (a : ℤ) : ℤ := rq a (2 ^ P)

/-- Fixed-point product. -/
def fmul (P : ℕ) (x y : ℤ) : ℤ := rs P (x * y)

/-- Fixed-point quotient (for `y > 0`). -/
def fdiv (P : ℕ) (x y : ℤ) : ℤ := rq (x * 2 ^ P) y

/-- Kronecker delta. -/
def kd (i k : ℕ) : ℝ := if i = k then 1 else 0

/-! ### Gaussian elimination without pivoting -/

/-- Fixed-point Gaussian elimination without pivoting, pivots in the order of `L`.
For `L = j :: L'`: multipliers `f_i = A_ij / A_jj`, Schur complement
`A'_ik = A_ik - f_i A_jk`, `b'_i = b_i - f_i b_j` on `L'`, recursive solve `z'`, and back
substitution `z_j = (b_j - Σ_{k∈L'} A_jk z'_k) / A_jj`.  Entries off `L` are `0`. -/
def geSolve (P : ℕ) : List ℕ → (ℕ → ℕ → ℤ) → (ℕ → ℤ) → ℕ → ℤ
  | [], _, _ => fun _ => 0
  | j :: L, A, b =>
    let f : ℕ → ℤ := fun i => fdiv P (A i j) (A j j)
    let z' := geSolve P L (fun i k => A i k - fmul P (f i) (A j k))
      (fun i => b i - fmul P (f i) (b j))
    fun k => if k = j then fdiv P (b j - ∑ k' ∈ L.toFinset, fmul P (A j k') (z' k')) (A j j)
      else z' k

section ge
variable (P : ℕ)

end ge

end J6Seg

end

section
-- Source module: Solutions.J6SegAlg
/-!
# Segmented inverse (S2): the algorithm as a pure function on fixed-point vectors

Everything here is integer arithmetic at precision `P` (an integer `a` stands for `a/2^P`);
the constants enter through a table `Tables` whose accuracy is a hypothesis of the error theorem
(`J6SegMain`).  Positions `0 ≤ i < s` are the indices of the resampling system `N = I + E`
(`Eexact`).  The blocks are runs `[st i, st i + len i)` of positions without wrap steps
(`Blocks`, validity `J6SegSplit.BlocksOK`); `w` is the coupling width.

Algorithm `segInv P T B U` (report L5-L9):
1. `Y = M⁻¹ U` block by block (`applyMinv`): seven fixed-point stages
   `D⁻¹, Γ⁻¹, L(a⁻), Γ, U(a⁺), Γ⁻¹, D` (`blkI`, q-LU of the recentred block).
2. Coupling of row `i`: offsets `0 < |h| ≤ w` leaving the block of `i` (`Hset`), column
   `col i h = (i + h) mod s`, entry `T.cpl i h ≈ 2^P e_{i,h}`.
3. Capacitance `C = I + Vᵀ M⁻¹ U` on the near-boundary positions `Jlist` (`capI`), using single
   entries of `M⁻¹` (`MinvI`, `TinvI` with the persymmetric flip so every entry is an `O(w)` sum).
4. `z = C⁻¹ (Vᵀ Y)` by elimination without pivoting in index order (`geSolve`).
5. `W = U z`, `Y₂ = M⁻¹ W`, output `X = Y - Y₂` (exact Woodbury).
-/

open Finset

namespace J6Seg

/-- Fixed-point constant tables (precision `P`). -/
structure Tables where
  /-- `≈ 2^P e^{-πα²θu²}` -/
  di : ℕ → ℤ
  /-- `≈ 2^P e^{πα²θu²}` (relative accuracy) -/
  dd : ℕ → ℤ
  /-- `≈ 2^P (q;q)_u` -/
  gam : ℕ → ℤ
  /-- `≈ 2^P / (q;q)_u` -/
  gi : ℕ → ℤ
  /-- `ap c h ≈ 2^P a⁺_h(β_c)` for a block starting at `c` -/
  ap : ℕ → ℕ → ℤ
  /-- `am c h ≈ 2^P a⁻_h(β_c)` -/
  am : ℕ → ℕ → ℤ
  /-- `cpl i h ≈ 2^P e_{i,h}` (coupling entries) -/
  cpl : ℕ → ℤ → ℤ

/-- Block structure: block of position `i` is `[st i, st i + len i)`; `w` = coupling width. -/
structure Blocks where
  s : ℕ
  st : ℕ → ℕ
  len : ℕ → ℕ
  w : ℕ

section alg
variable (P : ℕ)

/-- diagonal stage -/
def dgI (c x : ℕ → ℤ) : ℕ → ℤ := fun u => fmul P (c u) (x u)
/-- lower triangular Toeplitz stage (one rounding per output) -/
def lwI (c x : ℕ → ℤ) : ℕ → ℤ := fun u => rs P (∑ k ∈ range (u + 1), c (u - k) * x k)
/-- upper triangular Toeplitz stage (one rounding per output) -/
def upI (n : ℕ) (c x : ℕ → ℤ) : ℕ → ℤ := fun u => rs P (∑ k ∈ Ico u n, c (k - u) * x k)

/-- Block solve `M_c⁻¹ = D Γ⁻¹ U(a⁺) Γ L(a⁻) Γ⁻¹ D⁻¹` on a block starting at `a` of length `n`
(local indices `0 ≤ u < n`). -/
def blkI (T : Tables) (a n : ℕ) (z : ℕ → ℤ) : ℕ → ℤ :=
  dgI P T.dd (dgI P T.gi (upI P n (T.ap a) (dgI P T.gam (lwI P (T.am a)
    (dgI P T.gi (dgI P T.di z))))))

/-- `M⁻¹` (block diagonal) applied to a position vector. -/
def applyMinv (T : Tables) (B : Blocks) (x : ℕ → ℤ) : ℕ → ℤ :=
  fun i => blkI P T (B.st i) (B.len i) (fun u => x (B.st i + u)) (i - B.st i)

/-- Coupling offsets of row `i`: `0 < |h| ≤ w` with `i + h` outside the block of `i`. -/
def Hset (B : Blocks) (i : ℕ) : Finset ℤ :=
  ((range (2 * B.w + 1)).image (fun j : ℕ => (j : ℤ) - B.w)).filter
    (fun h => h ≠ 0 ∧ ¬ ((B.st i : ℤ) ≤ i + h ∧ (i : ℤ) + h < B.st i + B.len i))

/-- Column of the coupling entry `(i, h)`. -/
def col (s i : ℕ) (h : ℤ) : ℕ := (((i : ℤ) + h) % s).toNat

/-- `Σ_{max(u,v) ≤ k < n} a⁺_{k-u} Γ_k a⁻_{k-v}` in fixed point. -/
def TsumI (T : Tables) (a n u v : ℕ) : ℤ :=
  rs P (∑ k ∈ Ico (max u v) n, T.ap a (k - u) * fmul P (T.gam k) (T.am a (k - v)))

/-- Entry `(u, v)` of the block inverse `T⁻¹`; when both indices are far from the end the
persymmetric entry `(n-1-v, n-1-u)` is used, so the sum has at most `min(u,v) + 2w + 4` terms. -/
def TinvI (T : Tables) (w a n u v : ℕ) : ℤ :=
  if max u v + 2 * w + 3 < n then
    fmul P (fmul P (T.gi (n - 1 - v)) (T.gi (n - 1 - u))) (TsumI P T a n (n - 1 - v) (n - 1 - u))
  else fmul P (fmul P (T.gi u) (T.gi v)) (TsumI P T a n u v)

/-- Entry `(i, k)` of `M⁻¹` (same block): `D_u T⁻¹_{uv} D_v⁻¹`. -/
def MinvI (T : Tables) (B : Blocks) (i k : ℕ) : ℤ :=
  fmul P (fmul P (T.dd (i - B.st i)) (T.di (k - B.st i)))
    (TinvI P T B.w (B.st i) (B.len i) (i - B.st i) (k - B.st i))

/-- Near-boundary positions (the columns that carry coupling). -/
def nearB (B : Blocks) (j : ℕ) : Bool := decide (j < B.st j + B.w ∨ B.st j + B.len j < j + B.w + 1)

/-- The capacitance index list, in natural order. -/
def Jlist (B : Blocks) : List ℕ := (List.range B.s).filter (nearB B)

/-- Capacitance matrix `C = I + Vᵀ M⁻¹ U` (row `a`, column `b`, both in `Jlist`). -/
def capI (T : Tables) (B : Blocks) (a b : ℕ) : ℤ :=
  (if a = b then 2 ^ P else 0) + ∑ i ∈ Ico (B.st a) (B.st a + B.len a), ∑ h ∈ Hset B i,
    if col B.s i h = b then fmul P (MinvI P T B a i) (T.cpl i h) else 0

/-- `W = U z` (coupling applied to the capacitance solution). -/
def coupleI (T : Tables) (B : Blocks) (z : ℕ → ℤ) : ℕ → ℤ :=
  fun i => ∑ h ∈ Hset B i, fmul P (T.cpl i h) (z (col B.s i h))

/-- **The segmented inverse** `X ≈ 2^P N⁻¹ (U / 2^P)`. -/
def segInv (T : Tables) (B : Blocks) (U : ℕ → ℤ) : ℕ → ℤ :=
  let Y := applyMinv P T B U
  let z := geSolve P (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 (Gaussian-integer pair) version: real and imaginary parts separately
(`N` is a real matrix). -/
def segInvC (T : Tables) (B : Blocks) (U : ℕ → ℤ × ℤ) : ℕ → ℤ × ℤ :=
  fun i => (segInv P T B (fun j => (U j).1) i, segInv P T B (fun j => (U j).2) i)

end alg

end J6Seg

end

section
-- Source module: Solutions.J6SegToep
/-!
# Segmented inverse (S2): the recentred block and its explicit inverse

For a block of length `n` with recentring value `b` (`|b| ≤ 1/2`), `σ ≥ 1`, `θ = σ - 1`:
* `TF u k = e^{-πα²(σ(k-u)² + 2b(k-u))}` (the Toeplitz block, `= realToep`),
* `TinvF u v = Γ_u⁻¹ Γ_v⁻¹ Σ_{max(u,v) ≤ k < n} a⁺_{k-u} Γ_k a⁻_{k-v}` (its inverse, from the q-LU),
  with `Σ_k TF u k TinvF k v = δ_uv` (`TF_mul_TinvF`) and **persymmetry**
  `TinvF u v = TinvF (n-1-v) (n-1-u)` (`TinvF_persym`, report L6);
* `MF = D T D⁻¹`, `MinvF = D T⁻¹ D⁻¹` with `D = diag(e^{πα²θu²})` (report L4), `MF · MinvF = I`;
* the exact seven-stage pipeline `blkE` (scale by `D⁻¹`, `Γ⁻¹`, lower Toeplitz `a⁻`, `Γ`,
  upper Toeplitz `a⁺`, `Γ⁻¹`, `D`) equals `MinvF` applied to the input (`blkE_eq`).
* coefficient bounds: `|a^±_h| ≤ 1.0001 e^{-πα²θh}`, `Σ_h |a^±_h| ≤ 1.05`, `Γ_u ∈ (0,1]`,
  `Γ_u⁻¹ ≤ 1.0001`, `|TinvF u v| ≤ 1.1` (under `α ≥ 2`, `α²θ ≥ 2`).
-/

open Finset Matrix Real

namespace J6Seg
open J6QLU

/-! ### Real stage operators on `ℕ`-indexed vectors -/

/-- diagonal stage `u ↦ c_u x_u` -/
def dg (c x : ℕ → ℝ) : ℕ → ℝ := fun u => c u * x u
/-- lower triangular Toeplitz stage `u ↦ Σ_{k ≤ u} c_{u-k} x_k` -/
def lw (c x : ℕ → ℝ) : ℕ → ℝ := fun u => ∑ k ∈ range (u + 1), c (u - k) * x k
/-- upper triangular Toeplitz stage `u ↦ Σ_{u ≤ k < n} c_{k-u} x_k` -/
def up (n : ℕ) (c x : ℕ → ℝ) : ℕ → ℝ := fun u => ∑ k ∈ Ico u n, c (k - u) * x k

section consts
variable (α σ b : ℝ)

/-- `Γ_u = (q;q)_u` with `q = e^{-2πα²σ}` -/
noncomputable def gamR (u : ℕ) : ℝ := qpoch (qR α σ) u
/-- `Γ_u⁻¹` -/
noncomputable def giR (u : ℕ) : ℝ := (qpoch (qR α σ) u)⁻¹
/-- `D_u = e^{πα²θu²}`, `θ = σ - 1` -/
noncomputable def ddR (u : ℕ) : ℝ := exp (π * α ^ 2 * (σ - 1) * (u : ℝ) ^ 2)
/-- `D_u⁻¹ = e^{-πα²θu²}` -/
noncomputable def diR (u : ℕ) : ℝ := exp (-(π * α ^ 2 * (σ - 1) * (u : ℝ) ^ 2))

/-- the Toeplitz block entries -/
noncomputable def TF (u k : ℕ) : ℝ :=
  exp (-(π * α ^ 2 * (σ * ((k : ℝ) - u) ^ 2 + 2 * b * ((k : ℝ) - u))))

/-- the explicit inverse (q-LU) entries -/
noncomputable def TinvF (n u v : ℕ) : ℝ :=
  giR α σ u * giR α σ v *
    ∑ k ∈ Ico (max u v) n, aplus α σ b (k - u) * gamR α σ k * aminus α σ b (k - v)

/-- recentred block `D T D⁻¹` -/
noncomputable def MF (u k : ℕ) : ℝ := ddR α σ u * TF α σ b u k * diR α σ k
/-- its inverse `D T⁻¹ D⁻¹` -/
noncomputable def MinvF (n u k : ℕ) : ℝ := ddR α σ u * TinvF α σ b n u k * diR α σ k

/-- the exact block pipeline `D Γ⁻¹ U(a⁺) Γ L(a⁻) Γ⁻¹ D⁻¹` -/
noncomputable def blkE (n : ℕ) (v : ℕ → ℝ) : ℕ → ℝ :=
  dg (ddR α σ) (dg (giR α σ) (up n (aplus α σ b) (dg (gamR α σ) (lw (aminus α σ b)
    (dg (giR α σ) (dg (diR α σ) v))))))

end consts

/-! ### From matrices on `Fin n` to `ℕ`-indexed sums -/

section block
variable (α σ b : ℝ)

/-- `TinvM` = the q-LU inverse matrix. -/
noncomputable def TinvM (n : ℕ) : Matrix (Fin n) (Fin n) ℝ :=
  GamInv (qR α σ) n * upT (aplus α σ b) * Gam (qR α σ) n * lowT (aminus α σ b) * GamInv (qR α σ) n

end block

/-! ### Coefficient bounds -/

section bounds
variable {α σ b : ℝ}

end bounds

end J6Seg

end

section
-- Source module: Solutions.J6SegBlk
/-!
# Segmented inverse (S2): fixed-point error of one block solve and of single `M⁻¹` entries

* `sum_stage`: one rounded linear stage `rs P (Σ c̃_k x̃_k)` against the exact `Σ c_k x_k`.
* `blkI_err` (report L7/L10 block part): with tables accurate to `E` ulps (`dd` to relative
  `E` ulps), input `|z| ≤ 2^P B`, `Dm ≥ max_u D_u` and `8 K (B+1) ≤ 2^P`, `K = (n+1)E + 1`:
  `|blkI z / 2^P - M_c⁻¹ (z/2^P)| ≤ 9 Dm K (B+1) / 2^P`.
* `TinvI_err`, `MinvI_err`: single entries of `T⁻¹` (direct or persymmetric) and of `M⁻¹`.
-/

open Finset Real

namespace J6Seg
open J6QLU

section stage
variable (P : ℕ)

end stage

/-! ### The block solve -/

section blk
variable (P : ℕ) (α σ b : ℝ)

/-- Table accuracy on a block of length `n` starting at `a`, in ulps. -/
structure BlkTab (T : Tables) (a n : ℕ) (E : ℝ) : Prop where
  di : ∀ u < n, |(T.di u : ℝ) - 2 ^ P * diR α σ u| ≤ E
  dd : ∀ u < n, |(T.dd u : ℝ) - 2 ^ P * ddR α σ u| ≤ E * ddR α σ u
  gam : ∀ u < n, |(T.gam u : ℝ) - 2 ^ P * gamR α σ u| ≤ E
  gi : ∀ u < n, |(T.gi u : ℝ) - 2 ^ P * giR α σ u| ≤ E
  ap : ∀ h < n, |(T.ap a h : ℝ) - 2 ^ P * aplus α σ b h| ≤ E
  am : ∀ h < n, |(T.am a h : ℝ) - 2 ^ P * aminus α σ b h| ≤ E

variable {P α σ b}

end blk

end J6Seg

end

section
-- Source module: Solutions.J6SegRes
/-!
# Segmented inverse (S2): the exact Woodbury residual identity (abstract form)

Abstract setting on positions `ℕ` (a position set `I`): a block-diagonal matrix `M` given by its
blocks `bl i` and entries `Mx`, its block inverse `Mi`, and a sparse coupling `Δ` given row-wise
by offsets `Hs i`, columns `cl i h` and coefficients `cf i h`, with every coupling column in the
capacitance index set `J`.

`res_bound` (report L8 + L10, residual form): if the four fixed-point stages of the segmented
inverse are accurate,
* `y ≈ M⁻¹ u` (`ε₁`), `W ≈ Δ z` (`ε₃`), `y₂ ≈ M⁻¹ W` (`ε₂`),
* `z + (M⁻¹ Δ z)|_J ≈ y|_J` (`ε₄`, the capacitance solve),

then `x = y - y₂` satisfies `|((M + Δ) x - u)_i| ≤ A(ε₁+ε₂) + ε₃ + d(ε₄ + μ ε₃ + ε₂)`
(`A`, `d`, `μ` bound `‖M‖`, `‖Δ‖`, `‖M⁻¹‖` row-wise). No exact intermediate solution is needed.
`Miv_bound`: `‖M⁻¹‖ ≤ 1/(1-m)` from `M M⁻¹ = I` and `‖M - I‖ ≤ m < 1` (via `near_id`).
-/

open Finset

namespace J6Seg

section res
variable (bl : ℕ → Finset ℕ) (Mx : ℕ → ℕ → ℝ) (Hs : ℕ → Finset ℤ) (cl : ℕ → ℤ → ℕ)
  (cf : ℕ → ℤ → ℝ)

/-- `(M f)_i` for a block-diagonal `M`. -/
noncomputable def Mv (f : ℕ → ℝ) (i : ℕ) : ℝ := ∑ k ∈ bl i, Mx i k * f k

/-- `(Δ f)_i` for the sparse coupling. -/
noncomputable def Dv (f : ℕ → ℝ) (i : ℕ) : ℝ := ∑ h ∈ Hs i, cf i h * f (cl i h)

end res

end J6Seg

end

section
-- Source module: Solutions.KFMBaseExp
open Finset Complex

namespace IntMul.HvdH.ExpPf
open IntMul.HvdH

/-- A pair of integers as a Gaussian integer in `ℂ`. -/
def toC (p : ℤ × ℤ) : ℂ := (p.1 : ℂ) + (p.2 : ℂ) * I

end IntMul.HvdH.ExpPf

end

section
-- Source module: Solutions.KFMBaseRExp
open Complex

namespace IntMul.HvdH.RExpPf
open IntMul.HvdH

/-- A pair of integers as a Gaussian integer in `ℂ`. -/
def toC (p : ℤ × ℤ) : ℂ := (p.1 : ℂ) + (p.2 : ℂ) * I

end IntMul.HvdH.RExpPf

end

section
-- Source module: Solutions.KFMBaseEap
open Complex

namespace IntMul.HvdH.EPf
open IntMul.HvdH

/-- A pair of integers as a Gaussian integer in `ℂ`. -/
def toC (p : ℤ × ℤ) : ℂ := (p.1 : ℂ) + (p.2 : ℂ) * I

end IntMul.HvdH.EPf

end

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

namespace IntMul.HvdH.JPf
open IntMul.HvdH

/-- A pair of integers as a Gaussian integer in `ℂ`. -/
def toC (p : ℤ × ℤ) : ℂ := (p.1 : ℂ) + (p.2 : ℂ) * I

variable (s t α : ℕ)

/-- `x = α²θ = α²(t - s)/s`. -/
noncomputable def xθ : ℝ := ((α ^ 2 * (t - s) : ℕ) : ℝ) / s

section
variable {s t α}

end

end IntMul.HvdH.JPf

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

/-- Validity of a block structure for the resampling system with lengths `s < t`. -/
structure BlocksOK (s t : ℕ) (B : Blocks) : Prop where
  hs : B.s = s
  cover : ∀ i < s, B.st i ≤ i ∧ i < B.st i + B.len i ∧ B.st i + B.len i ≤ s
  cons : ∀ i < s, ∀ k, B.st i ≤ k → k < B.st i + B.len i → B.st k = B.st i ∧ B.len k = B.len i
  wf : ∀ i < s, ∀ j : ℕ, j + 1 < B.len i →
    bOf s t ((B.st i : ℤ) + j + 1) = bOf s t ((B.st i : ℤ) + j) + (t - s)
  hw : B.w ≤ s

/-- The block of position `i`. -/
def blS (B : Blocks) (i : ℕ) : Finset ℕ := Ico (B.st i) (B.st i + B.len i)

/-- Block part of `N`. -/
noncomputable def MxG (s t α : ℕ) (i k : ℕ) : ℝ := kd i k + eEx s t α i ((k : ℤ) - i)

/-- Exact block inverse `D T^{(b)}⁻¹ D⁻¹` (recentred, `σ = t/s`, `b = bOf(st)/s`). -/
noncomputable def MiG (s t α : ℕ) (B : Blocks) (i k : ℕ) : ℝ :=
  MinvF (α : ℝ) ((t : ℝ) / s) ((bOf s t (B.st i) : ℝ) / s) (B.len i) (i - B.st i) (k - B.st i)

/-- Coupling coefficients. -/
noncomputable def cfG (s t α : ℕ) (i : ℕ) (h : ℤ) : ℝ := eEx s t α i h

/-- `N x` on positions. -/
noncomputable def NxR (s t α : ℕ) (x : ℕ → ℝ) (i : ℕ) : ℝ :=
  x i + ∑' h : ℤ, eEx s t α i h * x ((((i : ℤ) + h) % s).toNat)

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

end glob

end J6Seg

end

section
-- Source module: Solutions.J6SegStage
/-!
# Segmented inverse (S2): accuracy of the global stages

With `K = (λ+1)E + 1` (`λ` ≥ every block length, `E` the table accuracy in ulps) and `Dm` ≥ every
`D_u = e^{πα²θu²}` used:
* `applyMinv_err`: `|applyMinv V / 2^P - M⁻¹ (V/2^P)| ≤ 9 Dm K (B+1)/2^P` (`|V/2^P| ≤ B`);
* `coupleI_err`: `|coupleI z / 2^P - Δ (z/2^P)| ≤ (2w+1)(1/2 + E Z)/2^P` (`|z/2^P| ≤ Z` on `J`);
* `capI_row`, `capI_apply`: the computed capacitance `C = I + (M⁻¹Δ)|_J` has row deviation
  `≤ 1/7 + λ(2w+1)δ` and `|C z - (z + M⁻¹ Δ z)| ≤ λ(2w+1) δ Z` with `δ = 15 Dm K / 2^P`.
-/

open Finset Real

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

/-- The hypotheses on the tables shared by all stages. -/
structure TabOK (P s t α : ℕ) (B : Blocks) (T : Tables) (E Dm : ℝ) : Prop where
  blk : ∀ i < s, BlkTab P (α : ℝ) ((t : ℝ) / s) ((bOf s t (B.st i) : ℝ) / s) T (B.st i) (B.len i) E
  cpl : ∀ i < s, ∀ h ∈ Hset B i, |(T.cpl i h : ℝ) - 2 ^ P * eEx s t α i h| ≤ E
  dm : ∀ i < s, ∀ u < B.len i, ddR (α : ℝ) ((t : ℝ) / s) u ≤ Dm

section stage
variable {P s t α : ℕ} {B : Blocks} {T : Tables} {E Dm : ℝ}

end stage

end J6Seg

end

section
-- Source module: Solutions.J6SegMain
/-!
# Segmented inverse (S2): the main error theorem (report L10)

`segInv_err`: for the resampling system `N = I + E` of lengths `s < t` with Gaussian parameter
`α ≥ 2` and `α²θ = α²(t-s)/s ≥ 3`, any valid block structure `B` (wrap-free runs of length
`≤ λ` partitioning `[0, s)`, coupling width `w ≤ s`), tables accurate to `E` ulps
(`TabOK`, `D_u ≤ Dm` on every block), a fixed-point input `U` (`|U/2^P| ≤ Bu`) approximating `v`
to `ε_in`, and the exact solution `w` of `N w = v`, the segmented inverse satisfies

`|segInv U i / 2^P - w_i| ≤ 60 G / 2^P + 12 e^{-πα² w(w+1)} (Bu+1) + (8/7) ε_in`

with `G = Dm ((λ+1)E+1) (Bu+1) (λ+1) (2w+1) (|J|+1)²` (`|J|` = number of near-boundary
positions), as soon as `1000 G ≤ 2^P`. In bits: precision
`P ≥ Q + log₂ Dm + 2 log₂(|J|+1) + log₂((λ+1)²(2w+1)) + log₂((λ+1)E+1) + log₂(Bu+1) + 16`
gives error `≤ 2^{-Q}` up to the truncation and input terms.
-/

open Finset Real

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

/-- The precision budget `G = Dm K (Bu+1) (λ+1) (2w+1) (|J|+1)²`. -/
noncomputable def segG (B : Blocks) (lam : ℕ) (E Dm Bu : ℝ) : ℝ :=
  Dm * (((lam : ℝ) + 1) * E + 1) * (Bu + 1) * ((lam : ℝ) + 1) * (2 * B.w + 1) *
    (((Jlist B).length : ℝ) + 1) ^ 2

section budget
variable {Dm K β L Wd Q e : ℝ}

end budget

end J6Seg

end