From c5612cbefef23f70f305542b99ce161d985c953b Mon Sep 17 00:00:00 2001
From: Ali Romman <226877337+aliromman@users.noreply.github.com>
Date: Mon, 26 Oct 2026 00:45:20 -0500
Subject: [PATCH 1/5] Retain exact finite-field arithmetic source for raw
 pencil extraction

---
 .../raw-pencil-decomposition/field_source.py  | 370 ++++++++++++++++++
 1 file changed, 370 insertions(+)
 create mode 100644 research/raw-pencil-decomposition/field_source.py

diff --git a/research/raw-pencil-decomposition/field_source.py b/research/raw-pencil-decomposition/field_source.py
new file mode 100644
index 0000000000..1066bb74c1
--- /dev/null
+++ b/research/raw-pencil-decomposition/field_source.py
@@ -0,0 +1,370 @@
+"""Exact finite-field tensor recognition. Ali Romman with OpenAI Codex. Apache-2.0."""
+from itertools import product
+from math import gcd, isqrt
+
+
+class Field:
+    """Use base-p digits and a monic polynomial, in constant-first order."""
+
+    def __init__(self, p, modulus):
+        if type(p) is not int or p < 2 or any(p % d == 0 for d in range(2, isqrt(p) + 1)):
+            raise ValueError('p must be prime')
+        if (not isinstance(modulus, (tuple, list)) or len(modulus) < 2
+                or modulus[-1] != 1
+                or any(type(c) is not int or not 0 <= c < p for c in modulus)):
+            raise ValueError('modulus must be monic with canonical prime-field coefficients')
+        self.p, self.e = p, len(modulus) - 1
+        self.modulus = tuple(modulus)
+        self.q = p ** self.e
+        digits = [[(a // p ** i) % p for i in range(self.e)] for a in range(self.q)]
+        self.add = [[sum(((x + y) % p) * p ** i for i, (x, y) in enumerate(zip(a, b)))
+                     for b in digits] for a in digits]
+        self.neg = [sum((-x % p) * p ** i for i, x in enumerate(a)) for a in digits]
+        self.mul = []
+        for a in digits:
+            row = []
+            for b in digits:
+                c = [0] * (2 * self.e - 1)
+                for i, x in enumerate(a):
+                    for j, y in enumerate(b):
+                        c[i + j] = (c[i + j] + x * y) % p
+                for k in range(len(c) - 1, self.e - 1, -1):
+                    for j in range(self.e):
+                        c[k - self.e + j] = (c[k - self.e + j] - c[k] * modulus[j]) % p
+                row.append(sum(c[i] * p ** i for i in range(self.e)))
+            self.mul.append(row)
+        self.inv = [0] + [self.power(a, self.q - 2) for a in range(1, self.q)]
+        if any(self.mul[a][self.inv[a]] != 1 for a in range(1, self.q)):
+            raise ValueError('modulus is reducible')
+
+    def power(self, a, n):
+        if n < 0:
+            raise ValueError('negative exponent')
+        r = 1
+        while n:
+            if n & 1:
+                r = self.mul[r][a]
+            a = self.mul[a][a]
+            n //= 2
+        return r
+
+    def total(self, xs):
+        r = 0
+        for x in xs:
+            r = self.add[r][x]
+        return r
+
+
+def identity(n):
+    return [[int(i == j) for j in range(n)] for i in range(n)]
+
+
+def columns(xs, n):
+    return [[v[i] for v in xs] for i in range(n)]
+
+
+def column(a, j):
+    return [r[j] for r in a]
+
+
+def rref(a, width, f):
+    rows = [r[:] for r in a]
+    pivots = []
+    k = 0
+    for j in range(width):
+        pivot = next((i for i in range(k, len(rows)) if rows[i][j]), None)
+        if pivot is None:
+            continue
+        rows[k], rows[pivot] = rows[pivot], rows[k]
+        rows[k] = [f.mul[f.inv[rows[k][j]]][x] for x in rows[k]]
+        for i in range(len(rows)):
+            if i != k and rows[i][j]:
+                c = f.neg[rows[i][j]]
+                rows[i] = [f.add[x][f.mul[c][y]] for x, y in zip(rows[i], rows[k])]
+        pivots.append(j)
+        k += 1
+        if k == len(rows):
+            break
+    return rows, pivots
+
+
+def rank(a, f):
+    return len(rref(a, len(a[0]) if a else 0, f)[1])
+
+
+def kernel(a, width, f):
+    rows, pivots = rref(a, width, f)
+    out = []
+    for j in range(width):
+        if j in pivots:
+            continue
+        v = [0] * width
+        v[j] = 1
+        for i, p in enumerate(pivots):
+            v[p] = f.neg[rows[i][j]]
+        out.append(v)
+    return out
+
+
+def inverse(a, f):
+    n = len(a)
+    rows, pivots = rref([r + e for r, e in zip(a, identity(n))], n, f)
+    if len(pivots) != n:
+        raise ValueError('singular matrix')
+    return [r[n:] for r in rows]
+
+
+def mv(a, v, f):
+    return [f.total(f.mul[x][y] for x, y in zip(r, v)) for r in a]
+
+
+def mm(a, b, f):
+    return columns([mv(a, column(b, j), f) for j in range(len(b[0]))], len(a))
+
+
+def determinant(a, f):
+    rows = [r[:] for r in a]
+    d = 1
+    for j in range(len(rows)):
+        k = next((i for i in range(j, len(rows)) if rows[i][j]), None)
+        if k is None:
+            return 0
+        if k != j:
+            rows[k], rows[j] = rows[j], rows[k]
+            d = f.neg[d]
+        c = rows[j][j]
+        d = f.mul[d][c]
+        for i in range(j + 1, len(rows)):
+            z = f.neg[f.mul[rows[i][j]][f.inv[c]]]
+            rows[i] = [f.add[x][f.mul[z][y]] for x, y in zip(rows[i], rows[j])]
+    return d
+
+
+def tensor_value(t, u, v, f):
+    n = len(t)
+    out = [0] * n
+    for i, x in enumerate(u):
+        if x:
+            for j, y in enumerate(v):
+                if y:
+                    c = f.mul[x][y]
+                    out = [f.add[z][f.mul[c][a]] for z, a in zip(out, t[i][j])]
+    return out
+
+
+def algebra_power(t, u, exponent, unit, f):
+    r = unit[:]
+    while exponent:
+        if exponent & 1:
+            r = tensor_value(t, r, u, f)
+        u = tensor_value(t, u, u, f)
+        exponent //= 2
+    return r
+
+
+def model(heights, f):
+    bounds = [f.p ** h for h in heights]
+    exponents = list(product(*(range(b) for b in bounds)))
+    ix = {v: i for i, v in enumerate(exponents)}
+    n = len(exponents)
+    t = [[[0] * n for _ in range(n)] for _ in range(n)]
+    for i, a in enumerate(exponents):
+        for j, b in enumerate(exponents):
+            c = tuple(x + y for x, y in zip(a, b))
+            if all(x < bound for x, bound in zip(c, bounds)):
+                t[i][j][ix[c]] = 1
+    return t, exponents
+
+
+def semilinear_kernel(t, unit, h, f):
+    """Return ker(a -> a^(p^h)), with the coefficient twist removed."""
+    n = len(t)
+    c = columns([algebra_power(t, v, f.p ** h, unit, f) for v in identity(n)], n)
+    untwist = f.p ** ((-h) % f.e)
+    return [[f.power(x, untwist) for x in v] for v in kernel(c, n, f)]
+
+
+def validate_tensor(t, f):
+    if not isinstance(f, Field):
+        raise ValueError('field contract required')
+    if not isinstance(t, (list, tuple)) or not t:
+        raise ValueError('positive square tensor required')
+    n = len(t)
+    if any(not isinstance(a, (list, tuple)) or len(a) != n for a in t):
+        raise ValueError('input dimensions differ')
+    for a in t:
+        for b in a:
+            if not isinstance(b, (list, tuple)) or len(b) != n:
+                raise ValueError('output dimension differs')
+            if any(type(x) is not int or not 0 <= x < f.q for x in b):
+                raise ValueError('noncanonical field entry')
+
+
+def check_certificate(t, cert, f):
+    """Check the factor maps by all basis products, without nucleus tests."""
+    n = len(t)
+    mu, _ = model(cert['heights'], f)
+    if len(mu) != n:
+        return False
+    u, v, w = (cert[k] for k in ('U', 'V', 'W'))
+    if any(determinant(a, f) == 0 for a in (u, v, w)):
+        return False
+    for i in range(n):
+        for j in range(n):
+            if tensor_value(t, column(u, i), column(v, j), f) != mv(w, mu[i][j], f):
+                return False
+    return True
+
+
+def recognize(t, f):
+    """Return (certificate, None), or (None, mathematical rejection reason)."""
+    validate_tensor(t, f)
+    n = len(t)
+    equations = []
+    for i, j, k in product(range(n), repeat=3):
+        row = [0] * (2 * n * n)
+        for a in range(n):
+            row[a * n + i] = f.add[row[a * n + i]][t[a][j][k]]
+            index = n * n + k * n + a
+            row[index] = f.add[row[index]][f.neg[t[i][j][a]]]
+        equations.append(row)
+    basis = kernel(equations, 2 * n * n, f)
+    if len(basis) != n:
+        return None, 'nucleus dimension'
+    xb = [[b[i * n:(i + 1) * n] for i in range(n)] for b in basis]
+    yb = [[b[n * n + i * n:n * n + (i + 1) * n] for i in range(n)] for b in basis]
+    if rank([b[:n*n] for b in basis], f) != n or rank([b[n*n:] for b in basis], f) != n:
+        return None, 'nucleus projection'
+    if any(mm(a, b, f) != mm(b, a, f) for a in yb for b in yb):
+        return None, 'noncommutative nucleus'
+    right = [[t[i][j][k] for j in range(n)] for i in range(n) for k in range(n)]
+    if rank(right, f) != n:
+        return None, 'right kernel'
+    flat = [sum(a, []) for a in yb]
+    _, positions = rref(flat, n * n, f)
+    coord_inverse = inverse([[a[j] for a in flat] for j in positions], f)
+
+    def coordinates(a):
+        v = sum(a, [])
+        c = mv(coord_inverse, [v[j] for j in positions], f)
+        if [f.total(f.mul[c[k]][flat[k][j]] for k in range(n)) for j in range(n*n)] != v:
+            raise ArithmeticError('nucleus closure')
+        return c
+
+    table = [[coordinates(mm(a, b, f)) for b in yb] for a in yb]
+    unit = coordinates(identity(n))
+    exponent = 1
+    while exponent < n:
+        exponent *= f.q
+    qmap = columns([algebra_power(table, v, exponent, unit, f) for v in identity(n)], n)
+    radical = kernel(qmap, n, f)
+    if len(radical) != n - 1:
+        return None, 'residue field or number of components'
+
+    def cyclic(module):
+        actions = [[[f.total(f.mul[c[k]][module[k][i][j]] for k in range(n))
+                     for j in range(n)] for i in range(n)] for c in radical]
+        jm = [column(a, j) for a in actions for j in range(n)]
+        if rank(jm, f) != n - 1:
+            return None
+        seed = next(v for v in identity(n) if rank(jm + [v], f) == n)
+        return columns([mv(a, seed, f) for a in module], n), seed
+
+    ux, wy = cyclic(xb), cyclic(yb)
+    if ux is None or wy is None:
+        return None, 'nonregular nucleus module'
+    u, u0 = ux
+    w, _ = wy
+    phi = columns([tensor_value(t, u0, e, f) for e in identity(n)], n)
+    v = mm(inverse(phi, f), w, f)
+    span = [tensor_value(table, a, b, f) for a in radical for b in radical]
+    current = rank(span, f)
+    generators, heights = [], []
+    h, power = 0, 1
+    while power < n:
+        h += 1
+        power *= f.p
+        for a in semilinear_kernel(table, unit, h, f):
+            new_rank = rank(span + [a], f)
+            if new_rank > current:
+                span.append(a)
+                generators.append(a)
+                heights.append(h)
+                current = new_rank
+    if f.p ** sum(heights) != n:
+        return None, 'height dimension'
+    mu, exponents = model(heights, f)
+    images = []
+    for exps in exponents:
+        value = unit[:]
+        for g, e in zip(generators, exps):
+            value = tensor_value(table, value, algebra_power(table, g, e, unit, f), f)
+        images.append(value)
+    c = columns(images, n)
+    if determinant(c, f) == 0:
+        raise ArithmeticError('height criterion did not produce a basis')
+    cert = {'heights': heights, 'U': mm(u, c, f), 'V': mm(v, c, f), 'W': mm(w, c, f)}
+    if not check_certificate(t, cert, f):
+        raise ArithmeticError('transported product mismatch')
+    index = gcd(f.q - 1, *(n * (f.p ** h - 1) // 2 for h in heights))
+    detu, detv, detw = (determinant(cert[k], f) for k in ('U', 'V', 'W'))
+    delta = f.mul[detw][f.inv[f.mul[detu][detv]]]
+    cert['sl_orbit_count'] = index
+    cert['sl_orbit_label'] = f.power(delta, (f.q - 1) // index)
+    return cert, None
+
+
+def bezout(a, b):
+    x, y, u, v = 1, 0, 0, 1
+    while b:
+        quotient, remainder = divmod(a, b)
+        a, b = b, remainder
+        x, u = u, x - quotient * u
+        y, v = v, y - quotient * v
+    return a, x, y
+
+
+def sl_transporter(first, second, f):
+    """Return determinant-one maps F,G,H with second(Fu,Gv)=H first(u,v)."""
+    a, why = recognize(first, f)
+    if a is None:
+        return None, 'first tensor: ' + why
+    b, why = recognize(second, f)
+    if b is None:
+        return None, 'second tensor: ' + why
+    if a['heights'] != b['heights'] or a['sl_orbit_label'] != b['sl_orbit_label']:
+        return None, 'different special-linear orbits'
+    n, m = len(first), f.q - 1
+    du, dv, dw = [f.mul[determinant(a[k], f)][f.inv[determinant(b[k], f)]]
+                  for k in ('U', 'V', 'W')]
+    z = f.mul[f.mul[du][dv]][f.inv[dw]]
+    primitive = next(c for c in range(1, f.q)
+                     if len({f.power(c, i) for i in range(m)}) == m)
+    logz = next(i for i in range(m) if f.power(primitive, i) == z)
+    g, coefficients = m, []
+    for h in a['heights']:
+        g, x, y = bezout(g, n * (f.p ** h - 1) // 2)
+        coefficients = [x * c for c in coefficients] + [y]
+    if logz % g:
+        raise ArithmeticError('orbit label disagrees with determinant image')
+    scalars = [f.power(primitive, (c * (logz // g)) % m) for c in coefficients]
+    _, exponents = model(a['heights'], f)
+    diagonal = [1] * n
+    for j, es in enumerate(exponents):
+        for c, e in zip(scalars, es):
+            diagonal[j] = f.mul[diagonal[j]][f.power(c, e)]
+    inverse_n = pow(n, -1, m)
+    unit_a = f.power(f.mul[du][f.inv[z]], inverse_n)
+    unit_b = f.power(f.mul[dv][f.inv[z]], inverse_n)
+    factors = [unit_a, unit_b, f.mul[unit_a][unit_b]]
+    out = {}
+    for key, target, scalar in zip(('U', 'V', 'W'), ('F', 'G', 'H'), factors):
+        middle = [[f.mul[scalar][diagonal[i]] if i == j else 0 for j in range(n)] for i in range(n)]
+        out[target] = mm(mm(b[key], middle, f), inverse(a[key], f), f)
+        if determinant(out[target], f) != 1:
+            raise ArithmeticError('determinant correction failed')
+    for i in range(n):
+        for j in range(n):
+            if tensor_value(second, column(out['F'], i), column(out['G'], j), f) != mv(out['H'], first[i][j], f):
+                raise ArithmeticError('special-linear product mismatch')
+    return out, None

From 2bd5bdb85039c65f0b5b314ed85153f8e6888a6f Mon Sep 17 00:00:00 2001
From: Ali Romman <226877337+aliromman@users.noreply.github.com>
Date: Mon, 26 Oct 2026 00:51:53 -0500
Subject: [PATCH 2/5] Extract complete raw finite-field pencil blocks and exact
 maps

---
 research/raw-pencil-decomposition/PROOF.md    | 165 ++++++++
 research/raw-pencil-decomposition/README.md   |  21 +
 research/raw-pencil-decomposition/RESULT.json |  56 +++
 .../raw-pencil-decomposition/SOURCES.json     |  25 ++
 research/raw-pencil-decomposition/TASK.json   |  64 ++++
 research/raw-pencil-decomposition/VERIFY.json |  20 +
 .../raw-pencil-decomposition/decompose.py     | 360 ++++++++++++++++++
 research/raw-pencil-decomposition/verify.py   | 196 ++++++++++
 8 files changed, 907 insertions(+)
 create mode 100644 research/raw-pencil-decomposition/PROOF.md
 create mode 100644 research/raw-pencil-decomposition/README.md
 create mode 100644 research/raw-pencil-decomposition/RESULT.json
 create mode 100644 research/raw-pencil-decomposition/SOURCES.json
 create mode 100644 research/raw-pencil-decomposition/TASK.json
 create mode 100644 research/raw-pencil-decomposition/VERIFY.json
 create mode 100644 research/raw-pencil-decomposition/decompose.py
 create mode 100644 research/raw-pencil-decomposition/verify.py

diff --git a/research/raw-pencil-decomposition/PROOF.md b/research/raw-pencil-decomposition/PROOF.md
new file mode 100644
index 0000000000..0034e0b2f7
--- /dev/null
+++ b/research/raw-pencil-decomposition/PROOF.md
@@ -0,0 +1,165 @@
+# Raw finite-field pencil decomposition
+
+Ali Romman, with substantial OpenAI Codex assistance. Apache-2.0.
+
+## Claim and interface
+
+Let K be a supplied finite field. For any two m by n matrices A and B over K,
+the algorithm returns the full Kronecker block list and invertible matrices
+R and Q such that A Q = R A0 and B Q = R B0. Thus the strict action is
+(R^-1, Q); the slice map is the identity. Zero rows, zero columns, zero-index
+blocks, dependent slices and the zero pencil are allowed. There is no need
+for an invertible slice over K.
+
+The output uses these models for s A0 + t B0:
+
+- Right index e: e by e+1, A0 has a unit diagonal and B0 a unit superdiagonal.
+- Left index e: the transpose of that right block.
+- Finite primary (f,e): (companion(f^e), I), with monic irreducible f.
+- Infinity of length e: (I, companion(x^e)).
+
+The code gives a deterministic algorithm polynomial in m+n for each fixed K.
+It uses field tables and exact row reduction. A coarse bound is
+O_K((m+n+1)^10) field operations and O_K((m+n+1)^4) live field entries.
+No bound polynomial in log(|K|), fast canonical-form bound, or numerical
+stability claim is made. The executable shares the usual finite memory and
+Python recursion limits. The written algorithm has no dimension bound.
+
+This completes the raw-input boundary for the supplied-model character
+interface. It does not compute the character formulas of that interface.
+It does not prove a physical supplier, native map, route, restoration schedule,
+precision bound, repeated preparation cost or complete kappa improvement.
+The selected .00310045 result and its main-transform bottleneck are unchanged.
+
+## 1. Shortest chains and singular blocks
+
+A polynomial right kernel vector of degree e has coefficients v0,...,ve with
+A v0=0, A vi+B v(i-1)=0 and B ve=0. Solve this homogeneous linear system in
+increasing e. If a right block exists, the first nonzero solution occurs at
+its least minimal index. Search through min(m,n-1); a right block of index e
+uses e rows and e+1 columns, so this range is complete.
+
+Use the classical Kronecker decomposition to justify the splitting step.
+A shortest chain has zero component on every regular or left block, and on
+every right block with larger index. There is no smaller right block. On each
+right block with the least index, it is a scalar multiple of its generator:
+a nonconstant multiplier would increase the degree. At least one scalar is
+nonzero. A change among equal-index copies makes this chain one direct
+summand. This also proves that its e+1 domain vectors and e image vectors
+are independent. No search through all vectors of K^n is required.
+
+Set ui=(-1)^i v(e-i), and wi=A ui for 0<=i<e. The chain equations give
+A ui=wi and B u(i+1)=wi, with the two endpoint zeros. These are precisely the
+stated right block. Let U and W contain these columns. Recover its complement
+by the linear equations
+
+```text
+F U = I, G W = I, G A = Ablock F, G B = Bblock F.
+```
+
+A solution exists because the embedding splits. The bases [U,ker F] and
+[W,ker G] put both coefficient matrices in block diagonal form. The code
+checks this identity before it recurses. For index zero, W has zero columns
+and G has zero rows; the same equations split an isolated zero column.
+
+After all right blocks are removed, apply the same step to the transpose.
+If At Qt=Rt Ct, then A Rt^-t=Qt^-t Ct^t. This gives the correct row and
+column maps for each left block, including an isolated zero row. The residue
+has no left or right polynomial kernel, hence is square and regular by
+Kronecker theory. This reasoning credits that theorem; it is not a new proof
+of its existence or uniqueness.
+
+## 2. Infinity without a rational invertible slice
+
+For the regular residue, start V0=0 and form V(i+1)={v:Bv is in A Vi}.
+Each step is a nullspace of [B,-A Vi] followed by projection onto the v
+coordinates. The increasing sequence stabilizes in at most n steps.
+
+On a finite block (C,I), every Vi is zero. On an infinity block (I,N),
+Vi=ker(N^i). Thus its limit is exactly the whole infinity space. A is
+injective there and its image W has the same dimension. B V is in W.
+Use W=A V to obtain the local pair (I,N), then solve the same linear
+projection equations to split this direct summand. The residual B is
+invertible. Normalize it to (B^-1 A,I). This proves the step even when each
+of the q+1 rational directions gives a singular matrix. No scalar search
+is used to select a global chart.
+
+## 3. Finite-field factorization and primary cyclic summands
+
+Find the minimal polynomial of T=B^-1 A from the first linear relation among
+I,T,...,T^n. Factor it with the classical Berlekamp method. For a squarefree
+polynomial f, the kernel of the q-power map minus the identity on K[x]/f
+has one dimension per irreducible factor. Its elements have values in K
+on each field factor. Successive gcds with b-c, for a basis b and all c in K,
+separate all factors. The code checks the final number of factors against
+that kernel dimension.
+
+For non-squarefree f, recurse on gcd(f,f') and its quotient. If f'=0,
+extract the p-th root of each coefficient with the inverse Frobenius power
+q/p and multiply the returned multiplicities by p. Each call reduces the
+degree. This includes extension fields and inseparable multiplicities.
+The verification uses a separate bounded monic-divisor enumeration, in
+addition to multiplying all factors back to the input polynomial.
+
+For each irreducible f with maximal exponent e, take ker(f(T)^e). These
+primary spaces have direct sum equal to the full space. On one primary space,
+find its nilpotency length e and a basis vector v with f(T)^(e-1)v nonzero.
+Its annihilator is exactly f^e. Its Krylov columns through degree deg(f)*e-1
+are independent and give the stated companion block. A cyclic submodule of
+maximal length in a direct sum of K[x]/f^j is a direct summand: a maximal
+component of v is a unit modulo f, so elementary module changes make it a
+generator in that component. Solve L U=I and L T=C L for its projection.
+Its kernel is a T-invariant complement. Repeat. Apply the same procedure
+with f=x to the infinity operator.
+
+All these steps are exact linear algebra or polynomial arithmetic with
+degrees bounded by m+n. The largest dense splitting system has O(N^2)
+unknowns and equations, with N=m+n+1. The stated coarse polynomial bound
+covers all degree searches, all splits, all primary steps and fixed-q
+Berlekamp splitting. It does not assert that this dense implementation is
+competitive with fast canonical-form algorithms.
+
+## 4. Certificate and independent acceptance boundary
+
+Both original equations A Q=R A0 and B Q=R B0 must hold and both maps must be
+invertible. Classification also requires irreducible primary labels and the
+complete block multiset; map equality alone would not justify a malformed
+primary label. The algorithm obtains those labels by the factorization above.
+The public `check` function checks maps to the requested model; it is not a
+standalone proof that arbitrary caller-supplied labels are canonical.
+
+The separate orbit oracle enumerates elementary invertible row and column
+operations on every raw pencil in five complete domains. Elementary additions
+and nonzero row/column scalings generate GL; swaps are products of these.
+It checks that the block type is constant on each orbit and differs between
+orbits. It does not use polynomial null chains, module decompositions or
+factorization to build its orbits. Every raw input also has both coefficient
+identities checked. Finite enumeration corroborates the algorithm; the
+all-dimension claim uses the arguments above and the classical module theorem.
+
+## Attribution and prior work
+
+[Nick Petalas PR #5250](https://github.com/CrocSwap/integer-mult-bounds/pull/5250)
+supplied the invertible-slice character formula and canonical-form reduction.
+[apoursh PR #5672](https://github.com/CrocSwap/integer-mult-bounds/pull/5672)
+supplied repeated singular blocks and explicit slice lifts.
+[PR #5683](https://github.com/CrocSwap/integer-mult-bounds/pull/5683) and
+[PR #5689](https://github.com/CrocSwap/integer-mult-bounds/pull/5689) supplied the
+mixed singular and complete finite-field model character interfaces. They
+explicitly leave raw decomposition separate. David Ganz's A030/A026/A029
+atlas lineage supplies the broader tensor-orbit motivation.
+
+The exact field and row-reduction module is retained from
+[PR #5664 source](https://github.com/aliromman/integer-mult-bounds/blob/8f7773e700d468396cbf893dfb18bf3ae6c20f66/research/residue-product-tensor-recognition/vendor/local.py).
+Its notice is unchanged. SOURCES.json pins its retained bytes. The remainder
+of this package is new integration, proof of this implementation and checks.
+
+Kronecker decomposition, primary module decomposition and Berlekamp
+factorization are classical. See [P. Van Dooren, 1979](https://perso.uclouvain.be/paul.vandooren/publications/VDooren79.pdf)
+for prior arbitrary-pencil extraction, and
+[Berlekamp's polynomial-factorization account](https://math.berkeley.edu/~berlek/poly.html).
+The preimage sequence is a Wong sequence; see
+[The Quasi-Kronecker Form for Matrix Pencils](https://epubs.siam.org/doi/abs/10.1137/110826278).
+We use their mathematical lineage, not their numerical implementations or
+complexity bounds. No global literature priority, human independent review,
+formal-kernel acceptance or new canonical-form theorem is claimed.
diff --git a/research/raw-pencil-decomposition/README.md b/research/raw-pencil-decomposition/README.md
new file mode 100644
index 0000000000..89b831389c
--- /dev/null
+++ b/research/raw-pencil-decomposition/README.md
@@ -0,0 +1,21 @@
+# Raw finite-field pencil decomposition
+
+This package extracts every singular and regular Kronecker block from two
+arbitrary rectangular finite-field matrices. It returns exact row and column
+maps and does not need an invertible base-field slice. It supplies the raw
+input boundary left by the supplied-model character packages.
+
+Read [PROOF.md](PROOF.md) for the algorithm, conventions, scope and credit.
+Read [RESULT.json](RESULT.json) for exact evidence and [SOURCES.json](SOURCES.json)
+for retained bytes. `decompose.decompose(A,B,m,n,field)` returns `(blocks,R,Q)`
+with `A Q=R A0` and `B Q=R B0`. Field elements use the retained module's
+base-p digit encoding. The slice map is the identity.
+
+```sh
+python3.12 -B scripts/research_package.py --package research/raw-pencil-decomposition
+```
+
+Expected: `PASS_RAW_PENCIL_DECOMPOSITION` and `PASS_REGISTERED_RESEARCH_PACKAGE`.
+Python 3.12 and the standard library suffice. No physical supplier or new
+kappa is claimed. The selected frontier is unchanged. Ali Romman with
+substantial OpenAI Codex assistance; no subagents. Apache-2.0.
diff --git a/research/raw-pencil-decomposition/RESULT.json b/research/raw-pencil-decomposition/RESULT.json
new file mode 100644
index 0000000000..f8072eb166
--- /dev/null
+++ b/research/raw-pencil-decomposition/RESULT.json
@@ -0,0 +1,56 @@
+{
+  "coefficient_equations": 22812,
+  "complete_domains": [
+    {
+      "pencils": 256,
+      "q": 2,
+      "shape": [
+        2,
+        2
+      ],
+      "strict_orbits": 16
+    },
+    {
+      "pencils": 6561,
+      "q": 3,
+      "shape": [
+        2,
+        2
+      ],
+      "strict_orbits": 24
+    },
+    {
+      "pencils": 65536,
+      "q": 4,
+      "shape": [
+        2,
+        2
+      ],
+      "strict_orbits": 34
+    },
+    {
+      "pencils": 4096,
+      "q": 2,
+      "shape": [
+        2,
+        3
+      ],
+      "strict_orbits": 20
+    },
+    {
+      "pencils": 4096,
+      "q": 2,
+      "shape": [
+        3,
+        2
+      ],
+      "strict_orbits": 20
+    }
+  ],
+  "constructed": 96,
+  "factor_polynomials": 1776,
+  "failure_controls": 183,
+  "no_invertible_slice": 12,
+  "source_inputs": 1,
+  "status": "PASS_RAW_PENCIL_DECOMPOSITION"
+}
diff --git a/research/raw-pencil-decomposition/SOURCES.json b/research/raw-pencil-decomposition/SOURCES.json
new file mode 100644
index 0000000000..69d2927c83
--- /dev/null
+++ b/research/raw-pencil-decomposition/SOURCES.json
@@ -0,0 +1,25 @@
+{
+  "schema_version": 1,
+  "base_commit": "df0137efd35e11b4445734bf612bdd4c7ebc330c",
+  "files": [
+    {
+      "path": "certificates/selected-result.json",
+      "usage": "provenance",
+      "git_blob": "581f8d87008354eb5b869f7ed9de1317ad837ce8",
+      "sha256": "2a0220e866417774606e0f991e41b9c1d26b1f43b97857238cf478dbc7d330b5"
+    },
+    {
+      "path": "docs/research/agent-task-catalog.json",
+      "usage": "provenance",
+      "git_blob": "f126dd2e2dad2b669e75d1d0d0bd05c600725602",
+      "sha256": "f7a48bfa2659b88378282a8eee7d9fe52cdf505594124c4d1686b58510fffd72"
+    },
+    {
+      "path": "research/raw-pencil-decomposition/field_source.py",
+      "usage": "input",
+      "git_blob": "1066bb74c19d180bfef7d9dd1c0650ec6bbd7780",
+      "sha256": "6daa160346d04691c6ca3c9c820b3df8982de4320ec01e5db10e3cb932f42c00",
+      "pin": "head:c5612cbefef23f70f305542b99ce161d985c953b:research/raw-pencil-decomposition/field_source.py"
+    }
+  ]
+}
diff --git a/research/raw-pencil-decomposition/TASK.json b/research/raw-pencil-decomposition/TASK.json
new file mode 100644
index 0000000000..f9104b8efa
--- /dev/null
+++ b/research/raw-pencil-decomposition/TASK.json
@@ -0,0 +1,64 @@
+{
+  "schema_version": 1,
+  "task_id": "bit-new-register-model",
+  "priority_revision": 7,
+  "base_commit": "df0137efd35e11b4445734bf612bdd4c7ebc330c",
+  "original_card": {
+    "id": "bit-new-register-model",
+    "question": "Can an explicitly larger or nonabelian register model escape the scoped three-register obstruction?",
+    "starting_points": [
+      "docs/research/research-family-status.md",
+      "research/bit-three-register-lower-bound/README.md",
+      "certificates/selected-result.json",
+      "research/bit-gl3-mixed-targets-length10/README.md",
+      "docs/research/community-efficient1000-review.md"
+    ],
+    "methods": [
+      "Construct the smallest exact identity in a newly specified legal operation alphabet.",
+      "Derive a target-specific relaxation bound for that same explicit alphabet."
+    ],
+    "acceptance_gate": "Name the changed hypothesis, all counted roles and a uniform identity; finite column tests alone do not certify it. A scoped local result may be retained without a complete supplier. Expanding it into further construction rounds requires the family-level budget and obligation gate in research-family-status.md; an optimistic passing relaxation alone is not acceptance.",
+    "first_step": {
+      "command": [
+        "python3",
+        "-B",
+        "scripts/research_baseline.py"
+      ],
+      "expected": "PASS_RESEARCH_BASELINE_BINDINGS_ONLY",
+      "inspect": "Read Model section 1 of the three-register bound and identify the exact excluded free operations.",
+      "scope": "Positive control / source-bound orientation only; this is not a new construction or full cold replay."
+    },
+    "first_experiment": "First compare the alphabet, role count and length with the fully replayed #1974 census and #2027 mixed-root tariff. Identify a changed premise and a legal native implementation; then price one smallest completed instance before searching larger words.",
+    "bottleneck_link": "A small-width completed saving changes the bit primitive; a finite-address identity alone is insufficient.",
+    "allowed_variations": [
+      "Choose a fixed larger register count or a specified nonabelian operation family."
+    ],
+    "scope_change_triggers": [
+      "Remaining inside the already-excluded three-register abelian model.",
+      "Repeating the two-register GL3 census or a covered three-register length without changing a hypothesis; a formal group-algebra word without native conversion is a partial result."
+    ],
+    "overlap_terms": [
+      "three-register",
+      "nonabelian",
+      "group algebra"
+    ]
+  },
+  "scope_status": "changed",
+  "actual_scope": "Complete raw finite-field rectangular matrix-pencil decomposition with explicit maps.",
+  "scope_change_reason": "The live tensor-atlas handoff leaves arbitrary-input extraction separate from the completed supplied-model characters.",
+  "bottleneck_link": "Provide canonical blocks and maps to evaluate two-slice tensor candidates; physical costs and the complete kappa bill remain open.",
+  "budget_minutes": 150,
+  "stop_rule": "Stop after full extraction proof and independent complete-domain map/type checks, or a counterexample.",
+  "outcome": "completed",
+  "evidence": "PROOF.md; RESULT.json; verify.py",
+  "author": "Ali Romman with substantial OpenAI Codex assistance",
+  "ai_assistance": "Substantial OpenAI Codex assistance; no subagents.",
+  "open_questions": [
+    {
+      "question": "Can a physical caller use an admissible raw pencil with a profitable complete bill?",
+      "why_it_matters": "Algebraic recognition alone does not improve kappa.",
+      "first_step": "Apply the supplied-model characters to these blocks and price every native map and route.",
+      "blocked_by": "No physical caller is supplied."
+    }
+  ]
+}
diff --git a/research/raw-pencil-decomposition/VERIFY.json b/research/raw-pencil-decomposition/VERIFY.json
new file mode 100644
index 0000000000..10e9abdd38
--- /dev/null
+++ b/research/raw-pencil-decomposition/VERIFY.json
@@ -0,0 +1,20 @@
+{
+  "schema_version": 1,
+  "status": "ready",
+  "python": "3.12",
+  "requirements": [],
+  "inputs": [
+    "research/raw-pencil-decomposition/field_source.py"
+  ],
+  "commands": [
+    {
+      "argv": [
+        "{python}",
+        "-B",
+        "{package}/verify.py"
+      ],
+      "timeout_seconds": 900,
+      "checks": "Complete strict GL orbits, literal maps, factor divisor oracle, mixed primaries, infinity and failure controls."
+    }
+  ]
+}
diff --git a/research/raw-pencil-decomposition/decompose.py b/research/raw-pencil-decomposition/decompose.py
new file mode 100644
index 0000000000..c8c275df66
--- /dev/null
+++ b/research/raw-pencil-decomposition/decompose.py
@@ -0,0 +1,360 @@
+"""Extract exact Kronecker blocks. Ali Romman with OpenAI Codex. Apache-2.0."""
+import field_source as a
+
+
+def tr(M, n):
+    return a.columns(M, n)
+
+
+def mm(A, B, n, f):
+    return a.columns([a.mv(A, a.column(B, j), f) for j in range(n)], len(A))
+
+
+def sub(A, r, c):
+    return [row[c:] for row in A[r:]]
+
+
+def span(vs, n, f):
+    rows, piv = a.rref(vs, n, f)
+    return rows[:len(piv)]
+
+
+def solve(rows, n, f):
+    rr, piv = a.rref(rows, n, f)
+    if any(not any(r[:n]) and r[n] for r in rr):
+        raise ValueError('inconsistent splitting equations')
+    x = [0] * n
+    for i, j in enumerate(piv):
+        x[j] = rr[i][n]
+    return x
+
+
+def split(A, B, m, n, U, W, C, D, b, c, f):
+    """Split an embedded summand of shape c by b with linear projections."""
+    size = b * n + c * m
+    eq = []
+    def add(terms, rhs):
+        row = [0] * (size + 1)
+        for j, x in terms:
+            row[j] = f.add[row[j]][x]
+        row[-1] = rhs
+        eq.append(row)
+    for i in range(b):
+        for j in range(b):
+            add([(i*n+k, U[k][j]) for k in range(n)], int(i == j))
+    for i in range(c):
+        for j in range(c):
+            add([(b*n+i*m+k, W[k][j]) for k in range(m)], int(i == j))
+    for X, Y in ((A, C), (B, D)):
+        for i in range(c):
+            for j in range(n):
+                add([(b*n+i*m+k, X[k][j]) for k in range(m)] +
+                    [(k*n+j, f.neg[Y[i][k]]) for k in range(b)], 0)
+    x = solve(eq, size, f)
+    F = [x[i*n:(i+1)*n] for i in range(b)]
+    G = [x[b*n+i*m:b*n+(i+1)*m] for i in range(c)]
+    Q = a.columns([a.column(U, j) for j in range(b)] + a.kernel(F, n, f), n)
+    R = a.columns([a.column(W, j) for j in range(c)] + a.kernel(G, m, f), m)
+    Ri = a.inverse(R, f)
+    X, Y = (mm(Ri, mm(Z, Q, n, f), n, f) for Z in (A, B))
+    for Z, model in ((X, C), (Y, D)):
+        if ([r[:b] for r in Z[:c]] != model or
+                any(any(r[b:]) for r in Z[:c]) or any(any(r[:b]) for r in Z[c:])):
+            raise ArithmeticError('split identity failed')
+    return X, Y, R, Q
+
+
+def chain(A, B, m, n, f):
+    for e in range(min(m, n - 1) + 1):
+        eq = []
+        for k in range(e + 2):
+            for i in range(m):
+                row = [0] * ((e + 1) * n)
+                if k <= e:
+                    row[k*n:(k+1)*n] = A[i]
+                if k:
+                    row[(k-1)*n:k*n] = B[i]
+                eq.append(row)
+        ker = a.kernel(eq, (e + 1) * n, f)
+        if ker:
+            v = ker[0]
+            U = a.columns([[f.mul[f.neg[1] if i % 2 else 1][x]
+                            for x in v[(e-i)*n:(e-i+1)*n]] for i in range(e+1)], n)
+            W = mm(A, U, e+1, f)
+            W = [r[:e] for r in W]
+            C = [[int(i == j) for j in range(e+1)] for i in range(e)]
+            D = [[int(i+1 == j) for j in range(e+1)] for i in range(e)]
+            return e, split(A, B, m, n, U, W, C, D, e+1, e, f)
+    return None
+
+
+def trim(p):
+    p = list(p)
+    while len(p) > 1 and p[-1] == 0:
+        p.pop()
+    return p
+
+
+def padd(p, q, f):
+    return trim([f.add[p[i] if i < len(p) else 0][q[i] if i < len(q) else 0]
+                 for i in range(max(len(p), len(q)))])
+
+
+def pmul(p, q, f):
+    r = [0] * (len(p)+len(q)-1)
+    for i, x in enumerate(p):
+        for j, y in enumerate(q):
+            r[i+j] = f.add[r[i+j]][f.mul[x][y]]
+    return trim(r)
+
+
+def pdiv(p, q, f):
+    p, q = trim(p), trim(q)
+    if q == [0]:
+        raise ValueError('zero divisor')
+    out = [0] * max(1, len(p)-len(q)+1)
+    while p != [0] and len(p) >= len(q):
+        j = len(p)-len(q)
+        c = f.mul[p[-1]][f.inv[q[-1]]]
+        out[j] = c
+        for i, x in enumerate(q):
+            p[i+j] = f.add[p[i+j]][f.neg[f.mul[c][x]]]
+        p = trim(p)
+    return trim(out), p
+
+
+def pgcd(p, q, f):
+    while q != [0]:
+        p, q = q, pdiv(p, q, f)[1]
+    return [f.mul[f.inv[p[-1]]][x] for x in p]
+
+
+def ppow(p, k, f, modulus=None):
+    r = [1]
+    while k:
+        if k & 1:
+            r = pmul(r, p, f)
+            if modulus is not None:
+                r = pdiv(r, modulus, f)[1]
+        p = pmul(p, p, f)
+        if modulus is not None:
+            p = pdiv(p, modulus, f)[1]
+        k //= 2
+    return r
+
+
+def factor(p, f):
+    """Use deterministic Berlekamp splitting, with inseparable recursion."""
+    p = trim(p)
+    if len(p) == 1:
+        return []
+    p = [f.mul[f.inv[p[-1]]][x] for x in p]
+    derivative = trim([f.mul[i % f.p][p[i]] for i in range(1, len(p))] or [0])
+    if derivative == [0]:
+        root = [f.power(p[i], f.q // f.p) for i in range(0, len(p), f.p)]
+        return sorted(factor(root, f) * f.p)
+    g = pgcd(p, derivative, f)
+    if len(g) > 1:
+        return sorted(factor(g, f) + factor(pdiv(p, g, f)[0], f))
+    d = len(p)-1
+    cols = []
+    for j in range(d):
+        v = ppow([0, 1], f.q*j, f, p)
+        v += [0] * (d-len(v))
+        v[j] = f.add[v[j]][f.neg[1]]
+        cols.append(v)
+    basis = a.kernel(a.columns(cols, d), d, f)
+    parts = [p]
+    for b in basis:
+        for c in range(f.q):
+            h = b[:]
+            h[0] = f.add[h[0]][f.neg[c]]
+            nxt = []
+            for part in parts:
+                g = pgcd(part, trim(h), f)
+                nxt += [g, pdiv(part, g, f)[0]] if 1 < len(g) < len(part) else [part]
+            parts = nxt
+    if len(parts) != len(basis):
+        raise ArithmeticError('Berlekamp split incomplete')
+    return sorted(tuple(x) for x in parts)
+
+
+def peval(p, T, f):
+    n = len(T)
+    out = [[0]*n for _ in range(n)]
+    for c in reversed(p):
+        out = mm(out, T, n, f)
+        for i in range(n):
+            out[i][i] = f.add[out[i][i]][c]
+    return out
+
+
+def minpoly(T, f):
+    n = len(T)
+    powers, P = [], a.identity(n)
+    for d in range(n+1):
+        flat = sum(P, [])
+        ker = a.kernel(a.columns(powers + [flat], n*n), d+1, f)
+        if ker:
+            v = next(x for x in ker if x[-1])
+            return [f.mul[f.inv[v[-1]]][x] for x in v]
+        powers.append(flat)
+        P = mm(P, T, n, f)
+    raise ArithmeticError('minimal polynomial absent')
+
+
+def companion(p, f):
+    n = len(p)-1
+    C = [[int(i == j+1) for j in range(n)] for i in range(n)]
+    for i in range(n):
+        C[i][-1] = f.neg[p[i]]
+    return C
+
+
+def primary(T, p, f):
+    n = len(T)
+    if not n:
+        return [], []
+    N = peval(p, T, f)
+    power, previous, e = a.identity(n), None, 0
+    while any(any(r) for r in power):
+        previous, power = power, mm(power, N, n, f)
+        e += 1
+    j = next(j for j in range(n) if any(a.column(previous, j)))
+    v = a.column(a.identity(n), j)
+    vs = []
+    d = (len(p)-1)*e
+    for _ in range(d):
+        vs.append(v)
+        v = a.mv(T, v, f)
+    U = a.columns(vs, n)
+    C = companion(ppow(p, e, f), f)
+    X, _, R, Q = split(T, a.identity(n), n, n, U, U, C, a.identity(d), d, d, f)
+    if R != Q:
+        raise ArithmeticError('operator split differs')
+    blocks, V = primary(sub(X, d, d), p, f)
+    return [(tuple(p), e)] + blocks, compose(Q, d, V, f)
+
+
+def operator(T, f):
+    n = len(T)
+    if not n:
+        return [], []
+    fac = factor(minpoly(T, f), f)
+    spaces, blocks = [], []
+    for p in sorted(set(fac)):
+        N = peval(ppow(p, fac.count(p), f), T, f)
+        basis = a.kernel(N, n, f)
+        V = a.columns(basis, n)
+        d = len(basis)
+        equations = []
+        for v in basis:
+            target = a.mv(T, v, f)
+            equations.append(solve([row + [x] for row, x in zip(V, target)], d, f))
+        restricted = a.columns(equations, d)
+        bs, U = primary(restricted, p, f)
+        spaces += [a.column(mm(V, U, d, f), j) for j in range(d)]
+        blocks += bs
+    return blocks, a.columns(spaces, n)
+
+
+def compose(Q, d, V, f):
+    n = len(Q)
+    E = [[int(i == j) if i < d and j < d else
+          V[i-d][j-d] if i >= d and j >= d else 0 for j in range(n)] for i in range(n)]
+    return mm(Q, E, n, f)
+
+
+def regular(A, B, f):
+    n = len(A)
+    basis = []
+    while True:
+        V = a.columns(basis, n)
+        AV = mm(A, V, len(basis), f)
+        eq = [r + [f.neg[x] for x in s] for r, s in zip(B, AV)]
+        nxt = span([v[:n] for v in a.kernel(eq, n+len(basis), f)], n, f)
+        if len(nxt) == len(basis):
+            break
+        basis = nxt
+    d = len(basis)
+    V = a.columns(basis, n)
+    W = mm(A, V, d, f)
+    BV = mm(B, V, d, f)
+    N = a.columns([solve([r+[x] for r, x in zip(W, a.column(BV, j))], d, f)
+                   for j in range(d)], d)
+    X, Y, R, Q = split(A, B, n, n, V, W, a.identity(d), N, d, d, f)
+    inf, U = primary(N, (0, 1), f)
+    residualB = sub(Y, d, d)
+    T = mm(a.inverse(residualB, f), sub(X, d, d), n-d, f)
+    finite, Z = operator(T, f)
+    localQ = [[U[i][j] if i < d and j < d else Z[i-d][j-d]
+               if i >= d and j >= d else 0 for j in range(n)] for i in range(n)]
+    BZ = mm(residualB, Z, n-d, f)
+    localR = [[U[i][j] if i < d and j < d else BZ[i-d][j-d]
+               if i >= d and j >= d else 0 for j in range(n)] for i in range(n)]
+    return [('infinity', (), e) for _, e in inf] + [('finite', p, e) for p, e in finite], mm(R, localR, n, f), mm(Q, localQ, n, f)
+
+
+def extract(A, B, m, n, f):
+    hit = chain(A, B, m, n, f)
+    if hit:
+        e, (X, Y, R, Q) = hit
+        c, b, tag = e, e+1, 'right'
+    else:
+        hit = chain(tr(A, n), tr(B, n), n, m, f)
+        if not hit:
+            if m != n:
+                raise ArithmeticError('regular residue is not square')
+            return regular(A, B, f)
+        e, (Xt, Yt, Rt, Qt) = hit
+        X, Y = tr(Xt, m), tr(Yt, m)
+        R, Q = tr(a.inverse(Qt, f), m), tr(a.inverse(Rt, f), n)
+        c, b, tag = e+1, e, 'left'
+    blocks, R2, Q2 = extract(sub(X, c, b), sub(Y, c, b), m-c, n-b, f)
+    return [(tag, (), e)] + blocks, compose(R, c, R2, f), compose(Q, b, Q2, f)
+
+
+def model(blocks, f):
+    pieces = []
+    for tag, p, e in blocks:
+        if tag in ('left', 'right'):
+            A = [[int(i == j) for j in range(e+1)] for i in range(e)]
+            B = [[int(i+1 == j) for j in range(e+1)] for i in range(e)]
+            m, n = e, e+1
+            if tag == 'left':
+                A, B, m, n = tr(A, n), tr(B, n), n, m
+        else:
+            C = companion(ppow(p if tag == 'finite' else (0, 1), e, f), f)
+            m = n = len(C)
+            A, B = (C, a.identity(n)) if tag == 'finite' else (a.identity(n), C)
+        pieces.append((A, B, m, n))
+    m, n = sum(x[2] for x in pieces), sum(x[3] for x in pieces)
+    A, B = [[0]*n for _ in range(m)], [[0]*n for _ in range(m)]
+    r = c = 0
+    for X, Y, h, w in pieces:
+        for i in range(h):
+            A[r+i][c:c+w], B[r+i][c:c+w] = X[i], Y[i]
+        r, c = r+h, c+w
+    return A, B, m, n
+
+
+def check(A, B, m, n, cert, f):
+    blocks, R, Q = cert
+    X, Y, h, w = model(blocks, f)
+    return ((h, w) == (m, n) and len(R) == m and all(len(r) == m for r in R)
+            and len(Q) == n and all(len(r) == n for r in Q)
+            and a.determinant(R, f) != 0 and a.determinant(Q, f) != 0
+            and mm(A, Q, n, f) == mm(R, X, n, f)
+            and mm(B, Q, n, f) == mm(R, Y, n, f))
+
+
+def decompose(A, B, m, n, f):
+    if not isinstance(f, a.Field) or any(type(x) is not int or x < 0 for x in (m, n)):
+        raise ValueError('explicit field and nonnegative dimensions required')
+    if any(len(M) != m or any(len(r) != n or any(type(x) is not int or not 0 <= x < f.q
+           for x in r) for r in M) for M in (A, B)):
+        raise ValueError('invalid coefficient matrix')
+    cert = extract(A, B, m, n, f)
+    if not check(A, B, m, n, cert, f):
+        raise ArithmeticError('coefficient identity failed')
+    return cert
diff --git a/research/raw-pencil-decomposition/verify.py b/research/raw-pencil-decomposition/verify.py
new file mode 100644
index 0000000000..afd693c62b
--- /dev/null
+++ b/research/raw-pencil-decomposition/verify.py
@@ -0,0 +1,196 @@
+#!/usr/bin/env python3
+"""Replay raw pencils against separate strict-equivalence orbits."""
+import copy
+import hashlib
+import itertools
+import json
+from pathlib import Path
+import random
+import sys
+
+ROOT = Path(__file__).resolve().parent
+sys.path.insert(0, str(ROOT))
+import decompose as d
+
+
+def need(ok, message):
+    if not ok:
+        raise ValueError(message)
+
+
+def authenticate():
+    sources = json.loads((ROOT/'SOURCES.json').read_text())
+    count = 0
+    for row in sources['files']:
+        if row['usage'] == 'input':
+            path = ROOT.parents[1]/row['path']
+            need(hashlib.sha256(path.read_bytes()).hexdigest() == row['sha256'], 'source hash')
+            count += 1
+    return count
+
+
+def flat_matrices(t, m, n):
+    return [list(t[i*n:(i+1)*n]) for i in range(m)], [list(t[m*n+i*n:m*n+(i+1)*n]) for i in range(m)]
+
+
+def neighbors(t, m, n, f):
+    for axis, size in ((0, m), (1, n)):
+        for i in range(size):
+            for j in range(size):
+                if i == j:
+                    continue
+                for c in range(1, f.q):
+                    out = list(t)
+                    for k in range(2):
+                        for z in range(n if axis == 0 else m):
+                            dst = k*m*n + (i*n+z if axis == 0 else z*n+i)
+                            src = k*m*n + (j*n+z if axis == 0 else z*n+j)
+                            out[dst] = f.add[out[dst]][f.mul[c][t[src]]]
+                    yield tuple(out)
+            for c in range(2, f.q):
+                out = list(t)
+                for k in range(2):
+                    for z in range(n if axis == 0 else m):
+                        dst = k*m*n + (i*n+z if axis == 0 else z*n+i)
+                        out[dst] = f.mul[c][out[dst]]
+                yield tuple(out)
+
+
+def literal(A, B, m, n, cert, f):
+    bs, R, Q = cert
+    X, Y, h, w = d.model(bs, f)
+    need((h, w) == (m, n), 'output shape')
+    need(d.a.determinant(R, f) != 0 and d.a.determinant(Q, f) != 0, 'basis rank')
+    for orig, canonical in ((A, X), (B, Y)):
+        for i, j in itertools.product(range(m), range(n)):
+            left = f.total(f.mul[orig[i][k]][Q[k][j]] for k in range(n))
+            right = f.total(f.mul[R[i][k]][canonical[k][j]] for k in range(m))
+            need(left == right, 'literal coefficient equation')
+
+
+def exhaustive(f, m, n):
+    seen, signatures = set(), set()
+    components = count = 0
+    for initial in itertools.product(range(f.q), repeat=2*m*n):
+        if initial in seen:
+            continue
+        todo = [initial]
+        seen.add(initial)
+        signature = None
+        while todo:
+            t = todo.pop()
+            A, B = flat_matrices(t, m, n)
+            cert = d.decompose(A, B, m, n, f)
+            literal(A, B, m, n, cert, f)
+            key = tuple(sorted(cert[0]))
+            if signature is None:
+                signature = key
+            need(key == signature, 'type changes inside independent GL orbit')
+            count += 1
+            for x in neighbors(t, m, n, f):
+                if x not in seen:
+                    seen.add(x)
+                    todo.append(x)
+        need(signature not in signatures, 'different GL orbits share extracted type')
+        signatures.add(signature)
+        components += 1
+    need(count == f.q**(2*m*n), 'complete raw domain')
+    row = {'q': f.q, 'shape': [m, n], 'pencils': count, 'strict_orbits': components}
+    print(json.dumps(row), flush=True)
+    return row
+
+
+def factor_oracle(f, degree):
+    cases = 0
+    for n in range(1, degree+1):
+        for coeff in itertools.product(range(f.q), repeat=n):
+            p = list(coeff) + [1]
+            factors = d.factor(p, f)
+            result = [1]
+            for g in factors:
+                result = d.pmul(result, g, f)
+                for k in range(1, (len(g)-1)//2+1):
+                    for cs in itertools.product(range(f.q), repeat=k):
+                        need(d.pdiv(g, list(cs)+[1], f)[1] != [0], 'factor reducible')
+            need(result == p, 'factor product')
+            cases += 1
+    return cases
+
+
+def random_basis(n, f, rng):
+    Q = d.a.identity(n)
+    for _ in range(5*n+1):
+        if n < 2:
+            break
+        i, j = rng.sample(range(n), 2)
+        c = rng.randrange(1, f.q)
+        Q[i] = [f.add[x][f.mul[c][y]] for x, y in zip(Q[i], Q[j])]
+    return Q
+
+
+def constructed():
+    rng = random.Random(7107)
+    count = equations = controls = no_slice = 0
+    for f in (d.a.Field(2, [0,1]), d.a.Field(3, [0,1]), d.a.Field(2, [1,1,1]), d.a.Field(3, [1,0,1])):
+        irreducibles = []
+        for k in (2, 3):
+            irreducibles.append(next(tuple(cs)+(1,) for cs in itertools.product(range(f.q), repeat=k)
+                                     if len(d.factor(list(cs)+[1], f)) == 1))
+        fixtures = [[], [('right',(),0)], [('left',(),0)],
+                    [('right',(),j) for j in range(5)], [('left',(),j) for j in range(5)],
+                    [('finite',irreducibles[0],3), ('finite',irreducibles[0],1),
+                     ('finite',irreducibles[1],2), ('infinity',(),3)],
+                    [('right',(),2),('right',(),2),('left',(),1),('left',(),3),
+                     ('finite',irreducibles[0],2),('finite',(0,1),2),('infinity',(),2)],
+                    [('finite',(f.neg[c],1),1) for c in range(f.q)] + [('infinity',(),1)]]
+        for bs in fixtures:
+            A, B, m, n = d.model(bs, f)
+            for _ in range(3):
+                R, Q = random_basis(m,f,rng), random_basis(n,f,rng)
+                X, Y = (d.mm(R,d.mm(Z,Q,n,f),n,f) for Z in (A,B))
+                cert = d.decompose(X,Y,m,n,f)
+                literal(X,Y,m,n,cert,f)
+                need(sorted(cert[0]) == sorted(bs), 'constructed block multiset')
+                count += 1
+                equations += 2*m*n
+                if bs == fixtures[-1]:
+                    for a, b in itertools.product(range(f.q), repeat=2):
+                        if a or b:
+                            Z = [[f.add[f.mul[a][x]][f.mul[b][y]] for x,y in zip(r,s)] for r,s in zip(X,Y)]
+                            need(d.a.determinant(Z,f) == 0, 'no invertible base-field slice')
+                    no_slice += 1
+                if m and n:
+                    bad = copy.deepcopy(cert)
+                    bad[1][0] = [0]*m
+                    need(not d.check(X,Y,m,n,bad,f), 'singular row-map control')
+                    bad = copy.deepcopy(cert)
+                    bad[2][0] = [0]*n
+                    need(not d.check(X,Y,m,n,bad,f), 'singular column-map control')
+                    changed = copy.deepcopy(X)
+                    changed[0][0] = f.add[changed[0][0]][1]
+                    need(not d.check(changed,Y,m,n,cert,f), 'coefficient corruption control')
+                    controls += 3
+    for A,B,m,n in [([[2]],[[0]],1,1), ([[0]],[[0,0]],1,1), ([[0]],[[0]],-1,1)]:
+        try:
+            d.decompose(A,B,m,n,d.a.Field(2,[0,1]))
+        except ValueError:
+            controls += 1
+        else:
+            raise ValueError('invalid-input control')
+    return {'constructed':count,'coefficient_equations':equations,'no_invertible_slice':no_slice,'failure_controls':controls}
+
+
+def main():
+    source_count = authenticate()
+    stats = [exhaustive(d.a.Field(p,mod),m,n) for p,mod,m,n in
+             [(2,[0,1],2,2),(3,[0,1],2,2),(2,[1,1,1],2,2),(2,[0,1],2,3),(2,[0,1],3,2)]]
+    factors = sum(factor_oracle(d.a.Field(p,mod),degree) for p,mod,degree in
+                  [(2,[0,1],7),(3,[0,1],5),(2,[1,1,1],4),(3,[1,0,1],3)])
+    result = {'status':'PASS_RAW_PENCIL_DECOMPOSITION','source_inputs':source_count,
+              'complete_domains':stats,'factor_polynomials':factors,**constructed()}
+    need(result == json.loads((ROOT/'RESULT.json').read_text()), 'committed result differs')
+    print(json.dumps(result,sort_keys=True))
+
+
+if __name__ == '__main__':
+    main()

From e9feb5299d7d22784c92a7b5db96a7aba4989db1 Mon Sep 17 00:00:00 2001
From: Ali Romman <226877337+aliromman@users.noreply.github.com>
Date: Mon, 26 Oct 2026 00:52:32 -0500
Subject: [PATCH 3/5] Include the registered research handoff section

---
 research/raw-pencil-decomposition/README.md | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/research/raw-pencil-decomposition/README.md b/research/raw-pencil-decomposition/README.md
index 89b831389c..e95fee4df0 100644
--- a/research/raw-pencil-decomposition/README.md
+++ b/research/raw-pencil-decomposition/README.md
@@ -19,3 +19,11 @@ Expected: `PASS_RAW_PENCIL_DECOMPOSITION` and `PASS_REGISTERED_RESEARCH_PACKAGE`
 Python 3.12 and the standard library suffice. No physical supplier or new
 kappa is claimed. The selected frontier is unchanged. Ali Romman with
 substantial OpenAI Codex assistance; no subagents. Apache-2.0.
+
+## Open next questions
+
+Connect the extracted blocks and maps to the supplied-model characters in
+PR5689, with an explicit convention check. For a physical candidate, price
+all maps, storage, preparation, routes and restoration before a kappa claim.
+The dense exact implementation is a reference algorithm; faster extraction
+can use established canonical-form methods with the same certificate boundary.

From 949b7de8d67dbd815cf12cc06d2758b6de4718da Mon Sep 17 00:00:00 2001
From: Ali Romman <226877337+aliromman@users.noreply.github.com>
Date: Mon, 26 Oct 2026 00:53:38 -0500
Subject: [PATCH 4/5] Use the registered replay timeout cap

---
 research/raw-pencil-decomposition/VERIFY.json | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/research/raw-pencil-decomposition/VERIFY.json b/research/raw-pencil-decomposition/VERIFY.json
index 10e9abdd38..70ad3d5b41 100644
--- a/research/raw-pencil-decomposition/VERIFY.json
+++ b/research/raw-pencil-decomposition/VERIFY.json
@@ -13,7 +13,7 @@
         "-B",
         "{package}/verify.py"
       ],
-      "timeout_seconds": 900,
+      "timeout_seconds": 600,
       "checks": "Complete strict GL orbits, literal maps, factor divisor oracle, mixed primaries, infinity and failure controls."
     }
   ]

From e34ae992744db135b9d4e378bfe26fd31fb82ad7 Mon Sep 17 00:00:00 2001
From: Ali Romman <226877337+aliromman@users.noreply.github.com>
Date: Mon, 26 Oct 2026 00:55:45 -0500
Subject: [PATCH 5/5] Record independent raw-pencil replay and portable
 evidence

---
 .../raw-pencil-decomposition/VALIDATION.md    | 32 +++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100644 research/raw-pencil-decomposition/VALIDATION.md

diff --git a/research/raw-pencil-decomposition/VALIDATION.md b/research/raw-pencil-decomposition/VALIDATION.md
new file mode 100644
index 0000000000..ecdf7d7ea8
--- /dev/null
+++ b/research/raw-pencil-decomposition/VALIDATION.md
@@ -0,0 +1,32 @@
+# Validation record
+
+Python3.14 normal replay and Python3.12.14 registered normal replay pass.
+A Python3.12.14 isolated optimized replay of the sealed science export also
+passes with no Git checkout and empty PATH. Ten exported files retain their
+exact bytes and read-only modes. An actual appended change to the retained
+field source fails authentication before mathematical checks.
+
+RESULT.json records 80545 raw pencils in five complete domains, 114 strict
+orbits, 1776 polynomial factorization cases, 96 constructed pencils, 22812
+constructed coefficient equations, 12 cases with no invertible rational slice
+and 183 failure controls. Every exhaustive raw pencil also has both maps
+checked. The independent GL orbit oracle covers all inputs in its stated
+domains. Larger dimensions are constructed evidence, not exhaustive evidence.
+
+The exact committed-byte and official export checks pass. The final export
+adds this record and the corrected 600-second manifest cap; the mathematical
+code and result are the sealed replay's bytes. Full-history ancestry preflight
+is omitted in this shallow checkout. Immediate source ancestry, the original
+remote field-source blob and all retained bytes were checked.
+
+The final overlap refresh covers all 7678 local records, 376 open PRs and
+5006 Discussion records with complete remote pagination, through
+2026-10-26T05:53:29Z. Canonical main remains
+`df0137efd35e11b4445734bf612bdd4c7ebc330c`. No competing raw decomposition
+claim was found. This does not cover unpublished work or every historical diff.
+
+The PR draft is checked with GitHub's Markdown renderer and its pinned file
+targets before publication. The published body and rendered HTML must then
+match that reviewed draft. Hosted CI, independent human review and proof-kernel
+acceptance remain separate. No selected-frontier edit, merge or closure of
+another contributor's work is part of this result.