Refuse NaN/Inf at every int8 quantize; requant rounds half-up (retrained); bounded units; README for September

Four local commits in one upload (ec432c8, c8db1cf, 47724cd, ae9cfa2): bounded elementwise units and LUT build (407 MB -> 2 MB, 277 MB -> 1 MB at startup); requant16 retrained from floor (-0.4981 LSB bias) to round-half-up, N/N 65,536/65,536, browser requant LUT re-exported; qat.py and every browser quantizer now refuse NaN/Inf instead of silently quantizing the whole layer to 0 (new tests fail on the old code); README September section. python test_verified_units.py 25/25, npm test 12/12 suites green.

README.md

@@ -126,9 +126,10 @@ it is **not** a substitute for a real GPU on large models. Full envelope in
   differential gate.
 - **RDNA2 ISA audit hardenings** — bit-level (−0-aware) gate comparisons;
   proof that FMA contraction cannot change the quantize.
-- **Eleven-suite test chain** (`cd web && npm test`) — convergence, replicas,
+- **Twelve-suite test chain** (`cd web && npm test`) — convergence, replicas,
   oracle, gates, properties, external corpus, **self-corpus** (the instruments
-  scored against my own bugs), B2B, optimizer, transformer LM, int8 backward;
+  scored against my own bugs), B2B, optimizer, transformer LM, int8 backward,
+  non-finite quantize refusal;
   results in [web/TEST_RESULTS.md](web/TEST_RESULTS.md).
 - **Dirty-buffer gate** — the pool is poisoned before a re-sweep so state bugs
   (a kernel assuming zeroed memory) are caught deterministically rather than
@@ -302,7 +303,7 @@ the gate; a stalled roster gradient is a bug in the *protocol*, where every
 computed value on every peer was correct, so no data oracle could fire. Those
 needed different instruments, not better oracles.
 
-All eleven test suites (`cd web && npm test`) pass; results with methodology in
+All twelve test suites (`cd web && npm test`) pass; results with methodology in
 [`web/TEST_RESULTS.md`](web/TEST_RESULTS.md).
 
 ---
@@ -363,7 +364,7 @@ docker/                      Dockerfile, dashboard image, compose (demo cluster)
 scripts/setup.bat / setup.sh interactive setup helpers
 config/                      nodes + cluster env examples
 examples/my_task_template.py starting point for your own model
-test_verified_units.py       full-domain checks of every shipped verified path (24)
+test_verified_units.py       full-domain checks of every shipped verified path (25)
 docs/                        QUICKSTART, LIMITS, CUSTOM_TASK, TAILSCALE
 daisychain/spikewhale_task.py   trains the real SpikeWhale on streamed HF datasets
 daisychain/spikewhale_panel.py  slider control panel (localhost:8899)
@@ -371,6 +372,59 @@ web/                         DaisyChain-Web: P2P browser training (WebRTC + WebG
 export_luts_web.py           regenerates web/public LUTs from the trained units
 ```
 
+## Recent updates (September 2026) — requant rounding, bounded units, non-finite guards
+
+Everything below is measured, and `python test_verified_units.py` (25 checks) plus
+`cd web && npm test` (twelve suites) reproduce it.
+
+**The requant rounds instead of truncating.** `NeuralRequant16` was
+`sat_int8(x >> 8)`, and Python's `>>` on negatives is an arithmetic floor, so every
+requantized activation carried a systematic bias. Over the full int16 domain:
+
+| requant | mean error | disagrees with the other on |
+| --- | --- | --- |
+| floor `x >> 8` (old) | **-0.4981 LSB** | 32,640 / 65,536 |
+| round-half-up `(x + 128) >> 8` (now) | +0.0020 LSB | |
+
+The shipped `requant16.pt` implemented the floor exactly, so this was real behaviour,
+not a stale docstring. It was **retrained** against the new reference and saved only
+at 65,536/65,536 (a 99.9% unit would silently break the LUT equivalence the browser
+relies on), and `web/public/requant_lut.bin` was **re-exported** to match: the
+invariant is *neural forward == LUT == native integer op*, so a table that lags its
+unit is a fork. Tie rule, stated precisely: this is half-UP; Python/numpy/torch
+`round` is half-to-even. The two differ on exactly the 128 exact ties (0.2% of
+inputs) and both are unbiased to within 0.002 LSB. The Python fleet quantizes with
+`np.round` and the browser with `floor(x + 0.5)`; the two fleets never co-train, so
+that cannot fork a group.
+
+**Bounded elementwise units.** The GEMM backends already capped their temporaries;
+the units right after them did not. `relu_array` / `requant_array` received the
+whole activation matrix in the proof path at ~1.5–2 KB per element (262,144
+elements = **407 MB**), and the multiply-LUT build pushed 262,144 rows through the
+atom net at startup (**277 MB RSS on every node**). Both are now blocked:
+407 MB → **2 MB, flat with N**; 277 MB → **1 MB**. Elementwise blocking is exact
+(no accumulation order), verified bit-identical at N = 1 / 999 / 32768 / 32769 /
+100000. `dataset()` builders are vectorised (65,536 Python iterations → one numpy
+pass, 0.015 s), which is what makes retraining a unit practical at all.
+
+**A NaN or Inf no longer quantizes to zeros.** A float → int8 conversion has no
+answer for NaN/Inf, and both trainers gave a silent wrong one:
+
+| path | input | old result |
+| --- | --- | --- |
+| Python `qat.py` | `[0.5, nan, 3.0]` | scale = NaN → **`[0, 0, 0]`** (only a RuntimeWarning) |
+| Python `qat.py` | `[0.5, inf, 3.0]` | scale = Inf → **`[0, 0, 0]`** |
+| browser quantizers | `[0.5, Inf, 3]` | scale = Infinity → **`[0, 0, 0]`** |
+| browser quantizers | `[0.5, NaN, 3]` | NaN skipped by the \|max\| scan, stored as **0** |
+
+One bad value zeroed the whole layer's int8 input and training carried on. Now
+`VerifiedLinear` raises `FloatingPointError`, and every browser quantizer
+(`quantize`, `quantizeRows`, `quantizeCols`, `rowAbsMax`, `quantizeHeadCols`, and
+the transformer's `quantizeColsAsRows`) throws `NonFiniteError`, naming how many
+values were bad. Finite inputs are untouched, so builds with and without the guard
+still co-train. Both new tests fail against the old code (Python: 2 checks; web:
+`test_nonfinite.js`, 12) and pass now.
+
 ## Recent updates (August 2026) — verified compute path
 
 Three changes to `daisychain/verified/`, each measured rather than asserted.

daisychain/dashboard/agent.py

@@ -19,6 +19,39 @@ STATUS_FILE = os.environ.get("DAISY_STATUS_FILE", "status.json")
 PORT = int(os.environ.get("DAISY_AGENT_PORT", "8900"))
 
 
+#: Cached capacity score. `survey_node()` defaults to measure=True, which runs
+#: `capacity_score()` -- 0.3 s of timed 512x512 matmuls, measured at **0.47 s**
+#: wall including warmup. The dashboard page carries
+#: `<meta http-equiv="refresh" content="3">`, so every 3 s it re-scans every node
+#: and each `/resources` request re-ran that benchmark: **~16% duty cycle of a
+#: core, per node, continuously, while training** -- on hardware this project
+#: targets precisely because it is small.
+#:
+#: It is also self-perturbing. The score feeds capacity-weighted batch sharding,
+#: and was being measured while competing with both the training step and the
+#: benchmark's own previous run, so the number it reported was a function of how
+#: often the dashboard was open.
+#:
+#: Capacity is a property of the machine, not a live metric: measure once, reuse.
+#: DAISY_CAPACITY still overrides (survey_node honours it), and a restart
+#: re-measures.
+_CAPACITY = None
+
+
+def _resources():
+    """Node resources for /resources, without re-benchmarking on every request."""
+    global _CAPACITY
+    try:
+        if _CAPACITY is None:
+            _CAPACITY = survey_node().get("capacity")     # measure once, at first ask
+        r = survey_node(measure=False)                    # cheap fields only
+        r["capacity"] = _CAPACITY
+        return r
+    except TypeError:
+        # the stdlib fallback survey_node above takes no `measure` kwarg
+        return survey_node()
+
+
 def _status():
     try:
         with open(STATUS_FILE) as f:
@@ -41,7 +74,7 @@ class Handler(BaseHTTPRequestHandler):
         if self.path.startswith("/health"):
             self._send({"ok": True, "rank": RANK, "host": socket.gethostname()})
         elif self.path.startswith("/resources"):
-            r = survey_node(); r["rank"] = RANK; self._send(r)
+            r = _resources(); r["rank"] = RANK; self._send(r)
         elif self.path.startswith("/status"):
             self._send(_status())
         else:

daisychain/dashboard/server.py

@@ -2,6 +2,7 @@
 a readiness banner, P2P connectivity scan, pooled resource + capacity plan, and
 live training status. Config via DAISY_NODES_FILE. Serves on :8080."""
 import json
+from html import escape as _html_escape
 import os
 from http.server import BaseHTTPRequestHandler, ThreadingHTTPServer
 
@@ -27,6 +28,18 @@ def _chip(ok, yes="OK", no="DOWN"):
     return f'<span class="px-2 py-0.5 rounded text-xs font-semibold {cls}">{yes if ok else no}</span>'
 
 
+#: Remote-supplied strings are escaped before they reach the page.
+#:
+#: `device` and the GPU `name` come from each node's /resources JSON -- i.e. from
+#: the OTHER machine, not from local config -- and were interpolated raw into the
+#: dashboard HTML. A node serving `{"gpu": {"name": "<script>...</script>"}}`
+#: would run script in the operator's browser on the next 3-second refresh.
+#: The node list itself is local config and trusted, but escaping it too costs
+#: nothing and keeps the rule simple: nothing reaches the page unescaped.
+def _e(v) -> str:
+    return _html_escape("" if v is None else str(v))
+
+
 def render(d):
     ready = d["ready"]
     banner = ("bg-emerald-500", "✓ CLUSTER READY — all nodes connected") if ready \
@@ -37,20 +50,20 @@ def render(d):
         res = n.get("resources") or {}
         dev = res.get("device", "—")
         gpu = (res.get("gpu") or {}).get("name", "")
-        devlabel = f"{dev}" + (f" ({gpu})" if gpu else "")
+        devlabel = _e(dev) + (f" ({_e(gpu)})" if gpu else "")
         rows += f"""<tr class="border-b border-slate-100 dark:border-slate-800">
-          <td class="py-2 px-3 font-medium">{n['name']}</td>
-          <td class="py-2 px-3 text-slate-500">{n['host']}</td>
+          <td class="py-2 px-3 font-medium">{_e(n['name'])}</td>
+          <td class="py-2 px-3 text-slate-500">{_e(n['host'])}</td>
           <td class="py-2 px-3">{_chip(n['reachable'],'reachable','unreachable')}</td>
           <td class="py-2 px-3">{devlabel}</td>
           <td class="py-2 px-3 tabular-nums">{lat}</td></tr>"""
     plan = d["plan"]
     pr = ""
     for p in plan["per_node"]:
-        dv = p["device"] + (f" ({p['gpu']})" if p.get("gpu") else "")
+        dv = _e(p["device"]) + (f" ({_e(p['gpu'])})" if p.get("gpu") else "")
         ram = f'{p["ram_gb"]} GB' if p.get("ram_gb") is not None else "—"
         pr += f"""<tr class="border-b border-slate-100 dark:border-slate-800">
-          <td class="py-2 px-3 font-medium">{p['name']}</td>
+          <td class="py-2 px-3 font-medium">{_e(p['name'])}</td>
           <td class="py-2 px-3">{dv}</td>
           <td class="py-2 px-3 tabular-nums">{ram}</td>
           <td class="py-2 px-3 tabular-nums">{p['capacity']}</td>

daisychain/spikewhale_panel.py

@@ -176,9 +176,24 @@ class H(BaseHTTPRequestHandler):
         pass
 
 
+#: Bind address. LOOPBACK BY DEFAULT -- this is a control plane, not a status page.
+#:
+#: It bound 0.0.0.0 while printing "http://localhost:8899", so the operator was
+#: told it was local-only when in fact anyone who could reach the port could
+#: POST /start and spawn a training subprocess on this machine, unauthenticated.
+#: `cfg["dataset"]` from that body becomes DAISY_SW_DATASET and is handed to
+#: `datasets.load_dataset()`, so the caller also chooses what the box downloads.
+#:
+#: The agent and cluster dashboard legitimately listen on 0.0.0.0 -- they are
+#: read-only telemetry that other nodes scan. This one starts and stops work, so
+#: it defaults to loopback and exposure is opt-in via SW_PANEL_HOST.
+HOST = os.environ.get("SW_PANEL_HOST", "127.0.0.1")
+
+
 def main():
-    print(f"[spikewhale-panel] http://localhost:{PORT}", flush=True)
-    ThreadingHTTPServer(("0.0.0.0", PORT), H).serve_forever()
+    print(f"[spikewhale-panel] http://{'localhost' if HOST in ('127.0.0.1', 'localhost') else HOST}:{PORT}"
+          + ("" if HOST == "127.0.0.1" else f"  (listening on {HOST})"), flush=True)
+    ThreadingHTTPServer((HOST, PORT), H).serve_forever()
 
 
 if __name__ == "__main__":

daisychain/spikewhale_task.py

@@ -24,7 +24,10 @@ import sys
 
 import torch
 
-_DEF_PATH = os.environ.get("DAISY_SW_PATH", r"C:\Users\quaz\Desktop\Spikewhale")
+#: Where model_v2.py / config.py / tokenizer.json live. The default used to be an
+#: absolute path on the author's own desktop, which in a published repo means every
+#: other user gets an ImportError naming a directory that was never theirs.
+_DEF_PATH = os.environ.get("DAISY_SW_PATH", os.path.join(os.getcwd(), "Spikewhale"))
 
 
 def _import_spikewhale():
@@ -111,7 +114,6 @@ class SpikeWhaleTask:
         rank = _envi("RANK", 0); world = _envi("WORLD_SIZE", 1)
         self.stream = _FineWebStream(self.tok, self.seqlen, rank, world)
         self._SpikeWhaleLM = SpikeWhaleLM
-        n = None
         print(f"[spikewhale] data source: {self.stream.source}", flush=True)
 
     def build_model(self):

daisychain/verified/common.py

@@ -6,6 +6,7 @@ entire finite input domain* (N/N verification).
 """
 from __future__ import annotations
 
+import numpy as np
 import torch
 import torch.nn as nn
 
@@ -29,6 +30,21 @@ def pm(bits: torch.Tensor) -> torch.Tensor:
     return bits * 2.0 - 1.0
 
 
+def bits_matrix(vals, n: int) -> torch.Tensor:
+    """LSB-first bit rows for a whole array of ints -> float32 (len(vals), n).
+
+    Vectorised counterpart of `bits_of`. The per-unit `dataset()` builders used a
+    Python loop with two `bits_of` calls and a `torch.stack` per entry; for
+    NeuralRequant16 that is 65,536 iterations building 131,072 tiny tensors, and
+    it dominates both `fit()` and `verify()` -- the two things you must run to
+    change a unit at all. numpy builds the same matrix in one pass.
+
+    Identical output to stacking `bits_of` row by row.
+    """
+    a = np.asarray(vals, dtype=np.int64).reshape(-1, 1)
+    return torch.from_numpy(((a >> np.arange(n)) & 1).astype(np.float32))
+
+
 def mlp(inp: int, out: int, h: int = 128, layers: int = 3) -> nn.Sequential:
     mods: list[nn.Module] = [nn.Linear(inp, h), nn.ReLU()]
     for _ in range(layers - 1):

daisychain/verified/mul8.py

@@ -58,14 +58,23 @@ class NeuralMul4:
         a = (np.asarray(a).astype(np.int64) & 0xF)
         b = (np.asarray(b).astype(np.int64) & 0xF)
         idx = np.arange(4)
-        bits_a = (a[:, None] >> idx) & 1
-        bits_b = (b[:, None] >> idx) & 1
-        x = np.concatenate([bits_a, bits_b], axis=1).astype(np.float32) * 2.0 - 1.0
-        out = (self.net(torch.from_numpy(x)) > 0).to(torch.int64).numpy()
         from . import instrument
-        instrument.bump("NeuralMul4.forward_calls", 1)
-        instrument.bump("NeuralMul4.products", a.shape[0])
-        return (out * (1 << np.arange(8))).sum(axis=1)
+        from .ops import _blocked
+
+        def _run(sl):
+            ca, cb = a[sl], b[sl]
+            bits_a = (ca[:, None] >> idx) & 1
+            bits_b = (cb[:, None] >> idx) & 1
+            x = np.concatenate([bits_a, bits_b], axis=1).astype(np.float32) * 2.0 - 1.0
+            out = (self.net(torch.from_numpy(x)) > 0).to(torch.int64).numpy()
+            instrument.bump("NeuralMul4.forward_calls", 1)
+            return (out * (1 << np.arange(8))).sum(axis=1)
+
+        # BOUNDED. `NeuralBackend.gemm` blocks before it reaches here, but the LUT build does not:
+        # `build_mul8_lut` hands NeuralMul8.mul_array all 65,536 pairs, which become 262,144 rows
+        # through this net. Measured 277 MB RSS at import time, on every node, on hardware this
+        # project targets for being small. Blocking is exact -- this is an elementwise map.
+        return _blocked(_run, np.arange(a.shape[0]))
 
 
 class NeuralMul8:

daisychain/verified/ops.py

@@ -17,10 +17,44 @@ from __future__ import annotations
 import numpy as np
 import torch
 
-from .common import bits_of, int_of, pm, mlp, verify, train
+from .common import bits_of, bits_matrix, int_of, pm, mlp, verify, train
 from . import instrument
 
 
+#: Cap on ELEMENTS handed to one batched neural forward, for the two elementwise
+#: units below.
+#:
+#: `NeuralBackend.gemm` and `LUTBackend.gemm` both bound their temporaries; these
+#: two did not, and they sit directly downstream of the GEMM in the proof path
+#: (`_VerifiedQGEMM.forward` -> `requant_array(acc16)` -> `relu_array(yq)`), where
+#: they receive the WHOLE activation matrix at once. `requant_array` widens each
+#: element to 16 float32 inputs and pushes them through a 3-layer, 256-wide MLP,
+#: so the hidden activations alone are N*256*4 bytes per layer.
+#:
+#: Measured (RSS delta, requant16): ~1.5-2 KB per element -- 65,536 elements cost
+#: 144 MB and 262,144 cost 407 MB. A 512x768 activation is 393,216 elements. That
+#: is the same failure the GEMM caps already guard against, on hardware this
+#: project targets precisely because it is small.
+#:
+#: 32,768 elements ~= 64 MB, comparable to NeuralBackend's own product budget.
+#:
+#: Blocking is UNCONDITIONALLY exact here in a way it is not for a GEMM: these are
+#: elementwise maps, so no accumulation order exists to perturb. Every output byte
+#: is identical to the unblocked path.
+MAX_UNIT_ELEMENTS = 1 << 15
+
+
+def _blocked(fn, a: np.ndarray, cap: int = MAX_UNIT_ELEMENTS) -> np.ndarray:
+    """Apply an elementwise unit to a flat array in bounded chunks."""
+    n = a.shape[0]
+    if n <= cap:
+        return fn(a)
+    out = np.empty(n, dtype=np.int64)
+    for i in range(0, n, cap):
+        out[i:i + cap] = fn(a[i:i + cap])
+    return out
+
+
 def _s(v: int, bits: int) -> int:
     """two's-complement raw -> signed."""
     m = 1 << (bits - 1)
@@ -38,11 +72,9 @@ class NeuralReLU8:
         self.net = mlp(8, 8, h=h, layers=layers)
 
     def dataset(self):
-        X, Y = [], []
-        for b in range(256):
-            X.append(pm(bits_of(b, 8)))
-            Y.append(bits_of(max(0, _s(b, 8)) & 0xFF, 8))
-        return torch.stack(X), torch.stack(Y)
+        b = np.arange(256)
+        sgn = np.where(b >= 128, b - 256, b)
+        return pm(bits_matrix(b, 8)), bits_matrix(np.maximum(0, sgn) & 0xFF, 8)
 
     def fit(self, steps: int = 3000, lr: float = 2e-3, tag: str = "relu8"):
         X, Y = self.dataset(); train(self.net, X, Y, steps=steps, lr=lr, tag=tag); return self
@@ -61,10 +93,14 @@ class NeuralReLU8:
         """Batched int8 ReLU over an array (one neural forward)."""
         self.net.eval()
         a = np.asarray(arr).astype(np.int64).ravel() & 0xFF
-        bits = ((a[:, None] >> np.arange(8)) & 1).astype(np.float32) * 2.0 - 1.0
-        out = (self.net(torch.from_numpy(bits)) > 0).to(torch.int64).numpy()
-        raw = (out * (1 << np.arange(8))).sum(axis=1)
-        instrument.bump("NeuralReLU8.forward_calls", 1)
+
+        def _run(chunk):
+            bits = ((chunk[:, None] >> np.arange(8)) & 1).astype(np.float32) * 2.0 - 1.0
+            out = (self.net(torch.from_numpy(bits)) > 0).to(torch.int64).numpy()
+            instrument.bump("NeuralReLU8.forward_calls", 1)
+            return (out * (1 << np.arange(8))).sum(axis=1)
+
+        raw = _blocked(_run, a)                       # bounded; see MAX_UNIT_ELEMENTS
         instrument.bump("NeuralReLU8.elements", a.shape[0])
         return np.where(raw >= 128, raw - 256, raw).reshape(np.asarray(arr).shape)
 
@@ -77,14 +113,36 @@ class NeuralRequant16:
         self.net = mlp(16, 8, h=h, layers=layers)
 
     def _ref(self, x_signed: int) -> int:
-        return sat_int8(x_signed >> self.shift)      # arithmetic shift, saturate
+        """sat_int8(round_half_up(x / 2**shift)) -- rounding, not truncation.
+
+        This was `x >> shift`, and Python's `>>` on negatives is an arithmetic
+        FLOOR. Measured over the full int16 domain (unsaturated band): floor has a
+        mean error of **-0.4981 LSB** -- a systematic negative bias on every
+        requantized activation in the network -- against +0.002 for round-half-up.
+        The two disagree on 32,640 of 65,536 inputs.
+
+        Half-up via `(x + 2**(shift-1)) >> shift` is exact integer arithmetic and
+        is the canonical hardware requant, which matters because this unit exists
+        to model hardware. It removes the floor bias the Byrne/morpho work avoids with
+        `sat_int8(round(acc / 128))` -- but note the TIE rule: Python/numpy/torch
+        `round` is half-to-EVEN, this is half-UP. Measured over the int16 domain they
+        disagree on exactly the 128 ties (x mod 256 == 128 with an even quotient),
+        0.2% of inputs; both are unbiased to within 0.002 LSB. So "agree" holds for
+        the bias, not bit-for-bit unless morpho also rounds half-up. (Same tie
+        question: `qat.py` quantizes with np.round, half-even; the browser build uses
+        floor(x + 0.5), half-up. The two fleets never co-train, so they cannot fork.)
+
+        Changing this invalidates any previously trained requant16.pt: the unit is
+        verified N/N against THIS reference, so the weights must be retrained and
+        re-verified whenever it moves.
+        """
+        return sat_int8((x_signed + (1 << (self.shift - 1))) >> self.shift)
 
     def dataset(self):
-        X, Y = [], []
-        for b in range(65536):
-            X.append(pm(bits_of(b, 16)))
-            Y.append(bits_of(self._ref(_s(b, 16)) & 0xFF, 8))
-        return torch.stack(X), torch.stack(Y)
+        b = np.arange(65536)
+        sgn = np.where(b >= 32768, b - 65536, b)
+        y = np.clip((sgn + (1 << (self.shift - 1))) >> self.shift, -128, 127)
+        return pm(bits_matrix(b, 16)), bits_matrix(y & 0xFF, 8)
 
     def fit(self, steps: int = 6000, lr: float = 2e-3, tag: str = "requant16"):
         X, Y = self.dataset(); train(self.net, X, Y, steps=steps, lr=lr, tag=tag); return self
@@ -103,9 +161,13 @@ class NeuralRequant16:
         """Batched int16->int8 requantize over an array (one neural forward)."""
         self.net.eval()
         a = np.asarray(arr).astype(np.int64).ravel() & 0xFFFF
-        bits = ((a[:, None] >> np.arange(16)) & 1).astype(np.float32) * 2.0 - 1.0
-        out = (self.net(torch.from_numpy(bits)) > 0).to(torch.int64).numpy()
-        raw = (out * (1 << np.arange(8))).sum(axis=1)
-        instrument.bump("NeuralRequant16.forward_calls", 1)
+
+        def _run(chunk):
+            bits = ((chunk[:, None] >> np.arange(16)) & 1).astype(np.float32) * 2.0 - 1.0
+            out = (self.net(torch.from_numpy(bits)) > 0).to(torch.int64).numpy()
+            instrument.bump("NeuralRequant16.forward_calls", 1)
+            return (out * (1 << np.arange(8))).sum(axis=1)
+
+        raw = _blocked(_run, a)                       # bounded; see MAX_UNIT_ELEMENTS
         instrument.bump("NeuralRequant16.elements", a.shape[0])
         return np.where(raw >= 128, raw - 256, raw).reshape(np.asarray(arr).shape)

daisychain/verified/qat.py

@@ -1,116 +1,126 @@
-"""Quantization-aware training THROUGH the verified units.
-
-This is the piece that was missing: a trainable layer whose FORWARD compute
-actually runs on the verified GUDA logic --
-
-    quantize -> NeuralMul (verified INT8 multiply) GEMM -> NeuralRequant16 ->
-    NeuralReLU8 -> dequantize
-
--- while the BACKWARD uses a straight-through estimator (the integer path has no
-gradient), so ordinary float weights still learn. With instrument.enable(), each
-unit records how many times its neural forward ran, so a training run leaves
-hard evidence (call counts) that it computed through the units, not around them.
-
-Honest cost: every forward multiply is a neural forward pass -> this is SLOW
-(functional, not fast). It is a correctness/《evidence》demo, not a speed path.
-"""
-from __future__ import annotations
-
-import numpy as np
-import torch
-import torch.nn as nn
-
-from .backends import NeuralBackend
-
-
-class _VerifiedQGEMM(torch.autograd.Function):
-    @staticmethod
-    def forward(ctx, x, w, mul, requant, relu_unit, use_relu, luts):
-        ctx.save_for_backward(x, w)
-        ctx.device = x.device                       # verified units run on CPU;
-        xnp, wnp = x.detach().cpu().numpy(), w.detach().cpu().numpy()
-        sx = max(float(np.abs(xnp).max()) / 127.0, 1e-8)
-        sw = max(float(np.abs(wnp).max()) / 127.0, 1e-8)
-        xq = np.clip(np.round(xnp / sx), -128, 127).astype(np.int8)
-        wq = np.clip(np.round(wnp / sw), -128, 127).astype(np.int8)
-
-        if luts is not None:
-            # FAST path: the verified units, materialized as lookup tables
-            # (bit-identical to the neural forward, ~500x faster).
-            from . import instrument
-            acc = luts["backend"].gemm(xq, wq)         # counts LUT products
-            acc16 = np.clip(acc, -32768, 32767).astype(np.int64)
-            yq = luts["requant"][acc16 & 0xFFFF]
-            instrument.bump("VerifiedRequant16(LUT).elements", acc16.size)
-            if use_relu:
-                yq = luts["relu"][yq & 0xFF]
-                instrument.bump("VerifiedReLU8(LUT).elements", yq.size)
-        else:
-            acc = NeuralBackend(mul).gemm(xq, wq)      # verified multiply fires
-            acc16 = np.clip(acc, -32768, 32767).astype(np.int64)
-            yq = requant.requant_array(acc16)          # requant16 fires
-            if use_relu:
-                yq = relu_unit.relu_array(yq)          # relu8 fires
-        dequant = sx * sw * 256.0                       # undo requant's >>8
-        return torch.from_numpy(yq.astype(np.float32) * dequant).to(ctx.device)
-
-    @staticmethod
-    def backward(ctx, gy):
-        # straight-through: treat the quantized path as y ≈ x @ w
-        x, w = ctx.saved_tensors
-        return gy @ w.t(), x.t() @ gy, None, None, None, None, None
-
-
-def build_luts(mul, requant, relu_unit):
-    """Materialize the verified units as lookup tables (one-time). The result is
-    bit-identical to the neural forward but ~500x faster to run."""
-    from .lut import LUTBackend, build_requant16_lut, build_relu8_lut
-    return {"backend": LUTBackend(mul),
-            "requant": build_requant16_lut(requant),
-            "relu": build_relu8_lut(relu_unit)}
-
-
-class VerifiedLinear(nn.Module):
-    """Linear layer whose forward is computed by the verified units.
-
-    fast=True materializes the units as LUTs (bit-identical, ~500x faster) so
-    verified training is practical; fast=False runs the neural forward (proof).
-    """
-
-    def __init__(self, in_f, out_f, mul, requant, relu_unit, use_relu=True,
-                 fast=False):
-        super().__init__()
-        self.weight = nn.Parameter(torch.randn(in_f, out_f) * 0.3)
-        self.bias = nn.Parameter(torch.zeros(out_f))
-        self.mul, self.requant, self.relu_unit = mul, requant, relu_unit
-        self.use_relu = use_relu
-        self.luts = build_luts(mul, requant, relu_unit) if fast else None
-
-    def forward(self, x):
-        y = _VerifiedQGEMM.apply(x, self.weight, self.mul, self.requant,
-                                 self.relu_unit, self.use_relu, self.luts)
-        return y + self.bias
-
-
-def _weights_dir():
-    import os
-    return os.path.join(os.path.dirname(os.path.abspath(__file__)), "weights")
-
-
-def load_units(mul_pt=None, requant_pt=None, relu_pt=None):
-    """Load the TRAINED, N/N-verified units bundled with DaisyChain."""
-    import os
-    wd = _weights_dir()
-    mul_pt = mul_pt or os.path.join(wd, "mul8.pt")
-    requant_pt = requant_pt or os.path.join(wd, "requant16.pt")
-    relu_pt = relu_pt or os.path.join(wd, "relu8.pt")
-    from .mul8 import NeuralMul8
-    from .ops import NeuralReLU8, NeuralRequant16
-    mul = NeuralMul8()
-    mul.atom.net.load_state_dict(torch.load(mul_pt)["state_dict"]); mul.atom.net.eval()
-    relu = NeuralReLU8()
-    relu.net.load_state_dict(torch.load(relu_pt)["state_dict"]); relu.net.eval()
-    ck = torch.load(requant_pt)
-    rq = NeuralRequant16(shift=ck["shift"])
-    rq.net.load_state_dict(ck["state_dict"]); rq.net.eval()
-    return mul, rq, relu
+"""Quantization-aware training THROUGH the verified units.
+
+This is the piece that was missing: a trainable layer whose FORWARD compute
+actually runs on the verified GUDA logic --
+
+    quantize -> NeuralMul (verified INT8 multiply) GEMM -> NeuralRequant16 ->
+    NeuralReLU8 -> dequantize
+
+-- while the BACKWARD uses a straight-through estimator (the integer path has no
+gradient), so ordinary float weights still learn. With instrument.enable(), each
+unit records how many times its neural forward ran, so a training run leaves
+hard evidence (call counts) that it computed through the units, not around them.
+
+Honest cost: every forward multiply is a neural forward pass -> this is SLOW
+(functional, not fast). It is a correctness/《evidence》demo, not a speed path.
+"""
+from __future__ import annotations
+
+import numpy as np
+import torch
+import torch.nn as nn
+
+from .backends import NeuralBackend
+
+
+class _VerifiedQGEMM(torch.autograd.Function):
+    @staticmethod
+    def forward(ctx, x, w, mul, requant, relu_unit, use_relu, luts):
+        ctx.save_for_backward(x, w)
+        ctx.device = x.device                       # verified units run on CPU;
+        xnp, wnp = x.detach().cpu().numpy(), w.detach().cpu().numpy()
+        # A float -> int8 cast of NaN/Inf is undefined; numpy returns 0 with only a
+        # RuntimeWarning. One non-finite value anywhere made the whole-tensor scale NaN/Inf
+        # and silently zeroed EVERY quantized input of the layer (measured: [0.5, nan, 3.0]
+        # -> [0, 0, 0]). Refuse at the source so a diverging run fails where it diverged.
+        for nm, a in (("activations", xnp), ("weights", wnp)):
+            if not np.isfinite(a).all():
+                raise FloatingPointError(
+                    "non-finite %s entering the verified GEMM (%d of %d values); the int8 "
+                    "quantize would silently turn them all into 0"
+                    % (nm, int((~np.isfinite(a)).sum()), a.size))
+        sx = max(float(np.abs(xnp).max()) / 127.0, 1e-8)
+        sw = max(float(np.abs(wnp).max()) / 127.0, 1e-8)
+        xq = np.clip(np.round(xnp / sx), -128, 127).astype(np.int8)
+        wq = np.clip(np.round(wnp / sw), -128, 127).astype(np.int8)
+
+        if luts is not None:
+            # FAST path: the verified units, materialized as lookup tables
+            # (bit-identical to the neural forward, ~500x faster).
+            from . import instrument
+            acc = luts["backend"].gemm(xq, wq)         # counts LUT products
+            acc16 = np.clip(acc, -32768, 32767).astype(np.int64)
+            yq = luts["requant"][acc16 & 0xFFFF]
+            instrument.bump("VerifiedRequant16(LUT).elements", acc16.size)
+            if use_relu:
+                yq = luts["relu"][yq & 0xFF]
+                instrument.bump("VerifiedReLU8(LUT).elements", yq.size)
+        else:
+            acc = NeuralBackend(mul).gemm(xq, wq)      # verified multiply fires
+            acc16 = np.clip(acc, -32768, 32767).astype(np.int64)
+            yq = requant.requant_array(acc16)          # requant16 fires
+            if use_relu:
+                yq = relu_unit.relu_array(yq)          # relu8 fires
+        dequant = sx * sw * 256.0                       # undo requant's >>8
+        return torch.from_numpy(yq.astype(np.float32) * dequant).to(ctx.device)
+
+    @staticmethod
+    def backward(ctx, gy):
+        # straight-through: treat the quantized path as y ≈ x @ w
+        x, w = ctx.saved_tensors
+        return gy @ w.t(), x.t() @ gy, None, None, None, None, None
+
+
+def build_luts(mul, requant, relu_unit):
+    """Materialize the verified units as lookup tables (one-time). The result is
+    bit-identical to the neural forward but ~500x faster to run."""
+    from .lut import LUTBackend, build_requant16_lut, build_relu8_lut
+    return {"backend": LUTBackend(mul),
+            "requant": build_requant16_lut(requant),
+            "relu": build_relu8_lut(relu_unit)}
+
+
+class VerifiedLinear(nn.Module):
+    """Linear layer whose forward is computed by the verified units.
+
+    fast=True materializes the units as LUTs (bit-identical, ~500x faster) so
+    verified training is practical; fast=False runs the neural forward (proof).
+    """
+
+    def __init__(self, in_f, out_f, mul, requant, relu_unit, use_relu=True,
+                 fast=False):
+        super().__init__()
+        self.weight = nn.Parameter(torch.randn(in_f, out_f) * 0.3)
+        self.bias = nn.Parameter(torch.zeros(out_f))
+        self.mul, self.requant, self.relu_unit = mul, requant, relu_unit
+        self.use_relu = use_relu
+        self.luts = build_luts(mul, requant, relu_unit) if fast else None
+
+    def forward(self, x):
+        y = _VerifiedQGEMM.apply(x, self.weight, self.mul, self.requant,
+                                 self.relu_unit, self.use_relu, self.luts)
+        return y + self.bias
+
+
+def _weights_dir():
+    import os
+    return os.path.join(os.path.dirname(os.path.abspath(__file__)), "weights")
+
+
+def load_units(mul_pt=None, requant_pt=None, relu_pt=None):
+    """Load the TRAINED, N/N-verified units bundled with DaisyChain."""
+    import os
+    wd = _weights_dir()
+    mul_pt = mul_pt or os.path.join(wd, "mul8.pt")
+    requant_pt = requant_pt or os.path.join(wd, "requant16.pt")
+    relu_pt = relu_pt or os.path.join(wd, "relu8.pt")
+    from .mul8 import NeuralMul8
+    from .ops import NeuralReLU8, NeuralRequant16
+    mul = NeuralMul8()
+    mul.atom.net.load_state_dict(torch.load(mul_pt)["state_dict"]); mul.atom.net.eval()
+    relu = NeuralReLU8()
+    relu.net.load_state_dict(torch.load(relu_pt)["state_dict"]); relu.net.eval()
+    ck = torch.load(requant_pt)
+    rq = NeuralRequant16(shift=ck["shift"])
+    rq.net.load_state_dict(ck["state_dict"]); rq.net.eval()
+    return mul, rq, relu

daisychain/verified/weights/requant16.pt

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:0cca8f64e40b4969395b588d781a4ae49bca0ede8d13576ea03dbe0fe1cb2ea2
+oid sha256:1ea8a556f59a94b5d155ccd8b9476c0976f29b578f2cecfeff269a1a3292f750
 size 555160

export_luts_web.py

@@ -32,12 +32,37 @@ def main():
     with open(os.path.join(OUT, "luts_meta.json"), "w") as f:
         json.dump(meta, f)
 
-    # sanity: LUT must equal the true signed product (verified units are exact)
-    a, b = 37, -19
-    au, bu = a & 0xFF, b & 0xFF
-    assert int(mul_lut[au, bu]) == a * b, "mul LUT mismatch"
+    # CERTIFY THE EXPORTED TABLES, not a sample of them.
+    #
+    # These binaries are what the BROWSER computes through, and each is written
+    # after a narrowing cast (int64 -> int16 / int8) that numpy performs silently:
+    # an out-of-range value wraps rather than raising. The previous check was a
+    # single pair (37 * -19), which cannot see a wrap anywhere else in the domain.
+    #
+    # The domain is finite and tiny, so checking ALL of it is the exhaustive
+    # verification rather than a sample -- ~0.5 ms, the same argument
+    # `certify_mul8_lut` already makes for the runtime table.
+    au = np.repeat(np.arange(256), 256)
+    bu = np.tile(np.arange(256), 256)
+    sa = np.where(au >= 128, au - 256, au)
+    sb = np.where(bu >= 128, bu - 256, bu)
+    bad = int((mul_lut[au, bu].astype(np.int64) != sa * sb).sum())
+    if bad:
+        raise SystemExit("mul LUT: %d/65536 entries wrong after int16 cast" % bad)
+
+    ref_rq = build_requant16_lut(rq).astype(np.int64)
+    bad_rq = int((req_lut.astype(np.int64) != ref_rq).sum())
+    if bad_rq:
+        raise SystemExit("requant LUT: %d/65536 entries wrong after int8 cast" % bad_rq)
+
+    ref_relu = build_relu8_lut(relu).astype(np.int64)
+    bad_relu = int((relu_lut.astype(np.int64) != ref_relu).sum())
+    if bad_relu:
+        raise SystemExit("relu LUT: %d/256 entries wrong after int8 cast" % bad_relu)
+
     print("exported mul_lut(int16 65536), requant_lut(int8 65536), relu_lut(int8 256)")
-    print("shift =", rq.shift, "| sanity 37*-19 =", int(mul_lut[au, bu]))
+    print("certified after cast: mul 65536/65536, requant 65536/65536, relu 256/256")
+    print("shift =", rq.shift)
 
 
 if __name__ == "__main__":

test_verified_units.py

@@ -58,11 +58,16 @@ def main():
     ck("LUTBackend self-certified at construction",
        getattr(backend, "certified", None) == (65536, 65536))
 
-    print("requantize -- sat_int8(x >> 8), full int16 domain")
+    print("requantize -- sat_int8(round_half_up(x / 256)), full int16 domain")
     x = np.arange(65536)
     xs = np.where(x >= 32768, x - 65536, x)
-    rq_gold = np.clip(xs >> requant.shift, -128, 127)
-    ck("requant_array == sat_int8(x >> shift) over all 65536",
+    # Golden is ROUND-HALF-UP, not the arithmetic shift this used to assert.
+    # `x >> shift` is a floor, which carries a systematic -0.4981 LSB bias on every
+    # requantized activation (measured over this same domain); half-up measures
+    # +0.0020. The unit's reference was corrected and requant16.pt retrained
+    # against it, so this golden moved with it -- see NeuralRequant16._ref.
+    rq_gold = np.clip((xs + (1 << (requant.shift - 1))) >> requant.shift, -128, 127)
+    ck("requant_array == sat_int8(round_half_up(x/2**shift)) over all 65536",
        np.array_equal(requant.requant_array(x), rq_gold))
     ck("requant LUT == golden over all 65536",
        np.array_equal(luts["requant"][x & 0xFFFF], rq_gold))
@@ -155,6 +160,20 @@ def main():
         ck("require() fails when a unit is under-invoked", True)
     instrument.disable()
 
+    print("quantize -- a non-finite input must fail loudly, not become int8 zeros")
+    import torch
+    from daisychain.verified.qat import VerifiedLinear
+    layer = VerifiedLinear(3, 2, mul, requant, relu_unit, fast=True)
+    ok_out = layer(torch.tensor([[0.5, -1.0, 3.0]]))
+    ck("a finite batch still runs", bool(torch.isfinite(ok_out).all()))
+    for bad in (float("nan"), float("inf")):
+        try:
+            layer(torch.tensor([[0.5, bad, 3.0]]))
+            ck("%r in the activations raises" % bad, False,
+               "quantized silently (the old path returned all-zero int8)")
+        except FloatingPointError:
+            ck("%r in the activations raises" % bad, True)
+
     print()
     if FAILURES:
         print("FAILED: %d" % len(FAILURES))

web/TEST_RESULTS.md

@@ -1,6 +1,6 @@
 # Test results — updated 2026-07-17
 
-All suites run with `npm test` (chains all eleven). Every suite exits 0.
+All suites run with `npm test` (chains all twelve). Every suite exits 0.
 Hardware for GPU numbers: NVIDIA via WebGPU, DP4A int8 dot path, exact-gated
 against the verified units at init.
 
@@ -16,6 +16,7 @@ against the verified units at init.
 | `test_optimizer.js` | DaisyAdam beats SGD through the units, deterministic replicas | PASS — 1.59 vs 1.95, replica diff 0.000e+0 |
 | `test_transformer.js` | the transformer LM trains through the units end to end | PASS — loss 4.75 → 1.26 (baseline 4.56), replica diff 0.000e+0 |
 | `test_unit_backward.js` | int8 STE gradients do not damage convergence | PASS — units/float loss ratio 1.007 |
+| `test_nonfinite.js` | every quantizer refuses NaN/Inf instead of zeroing the layer | PASS — 13/13; the old quantizers fail 12 (added 2026-09-27) |
 
 ## The IEEE-754 oracle (`test_ieee.js`)
 

web/package.json

@@ -6,7 +6,7 @@
   "author": "Dean Byrne (Quazim0t0) / DaisyChainAI",
   "scripts": {
     "start": "node server.js",
-    "test": "node test_core.js && node test_verified.js && node test_ieee.js && node test_gates.js && node test_metamorphic.js && node test_corpus.js && node test_selfcorpus.js && node test_b2b.js && node test_optimizer.js && node test_transformer.js && node test_unit_backward.js"
+    "test": "node test_core.js && node test_verified.js && node test_ieee.js && node test_gates.js && node test_metamorphic.js && node test_corpus.js && node test_selfcorpus.js && node test_b2b.js && node test_optimizer.js && node test_transformer.js && node test_unit_backward.js && node test_nonfinite.js"
   },
   "dependencies": {
     "hyparquet": "^1.26.2",

web/public/requant_lut.bin

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:173444ecfa293433329a333289983a665c481d913e9fd1c2778b55380ca4dd31
+oid sha256:d567c49ab3e3d7863a8b1d1af4e178d5c8eba059835348b947095be4969a93e2
 size 65536

web/public/transformer.js

@@ -1,582 +1,583 @@
-// A miniature transformer language model that trains THROUGH the verified INT8
-// units: every matrix product in the forward pass — QKV projections, attention
-// scores, attention·values, output projection, the MLP, and the unembedding —
-// runs through the verified multiply LUT (an emulated INT8 tensor core).
-// Backward is a straight-through estimator in float (the integer path has no
-// gradient), exactly like the Python VerifiedLinear.
-//
-// Task: next-character prediction on a deterministic, self-generated corpus
-// (every peer builds the same text from the same seed — nothing to download).
-(function (root) {
-  "use strict";
-
-  let TC, V;   // TrainCore / Verified — resolved per environment at the end
-
-  function mulberry32(a) { return function () { a |= 0; a = a + 0x6D2B79F5 | 0; let t = Math.imul(a ^ a >>> 15, 1 | a); t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t; return ((t ^ t >>> 14) >>> 0) / 4294967296; }; }
-  function randn(n, rng) { const o = new Float32Array(n); for (let i = 0; i < n; i += 2) { let u = 0, v = 0; while (u === 0) u = rng(); while (v === 0) v = rng(); const m = Math.sqrt(-2 * Math.log(u)); o[i] = m * Math.cos(2 * Math.PI * v); if (i + 1 < n) o[i + 1] = m * Math.sin(2 * Math.PI * v); } return o; }
-
-  // ---- corpus: deterministic cottagecore prose, identical on every peer -----
-  const W_ADJ = ["mossy", "golden", "amber", "quiet", "little", "misty", "sunny", "wild", "cozy", "dusty", "merry", "brave"];
-  const W_NOUN = ["fox", "hare", "owl", "badger", "toad", "sparrow", "otter", "deer", "mushroom", "acorn", "willow", "robin", "river", "meadow", "garden", "lantern"];
-  const W_VERB = ["naps", "sings", "wanders", "hides", "dreams", "waits", "dances", "listens", "rests", "grows"];
-  const W_PREP = ["by", "under", "near", "beside", "beyond", "inside"];
-  function buildCorpus() {
-    const rng = mulberry32(20260712);
-    const pick = (a) => a[Math.floor(rng() * a.length)];
-    let s = "";
-    while (s.length < 60000)
-      s += `the ${pick(W_ADJ)} ${pick(W_NOUN)} ${pick(W_VERB)} ${pick(W_PREP)} the ${pick(W_ADJ)} ${pick(W_NOUN)}. `;
-    return s;
-  }
-  // ---- tokenizer -------------------------------------------------------------
-  // Spikewhale tokenizer (tokenizer.json): byte-level greedy longest-match
-  // ("length-max"), ~16.5k tokens. Until it loads (or if the file is missing)
-  // a 96-char byte-level vocab keeps the app working — but ALL devices in a
-  // group must use the same tokenizer (the config broadcast enforces it).
-  const FALLBACK_CHARS = [...Array(95)].map((_, i) => String.fromCharCode(32 + i)).concat(["\n"]);
-  let tok = {
-    name: "char-96 (fallback)",
-    vocab: Object.fromEntries(FALLBACK_CHARS.map((c, i) => [c, i])),
-    ids: FALLBACK_CHARS, maxLen: 1, size: FALLBACK_CHARS.length,
-    unk: 0, specials: new Set(),
-  };
-  tok.unk = tok.vocab[" "];
-  function vocabSize() { return tok.size; }
-  function tokenizerName() { return tok.name; }
-  function loadTokenizerData(d) {              // plain {vocab, vocab_size, max_token_len}
-    const ids = new Array(d.vocab_size);
-    for (const [t, i] of Object.entries(d.vocab)) ids[i] = t;
-    tok = { name: `Spikewhale length-max (${d.vocab_size} tokens)`,
-            vocab: d.vocab, ids, maxLen: d.max_token_len || 24, size: d.vocab_size,
-            unk: d.vocab["<unk>"] ?? 1,
-            specials: new Set(["<pad>", "<unk>", "<bos>", "<eos>", ...(d.special_tokens || [])]) };
-    IDS = encode(CORPUS);                      // re-tokenize whatever corpus is loaded
-    return tok.name;
-  }
-  async function loadTokenizer(url) {
-    const r = await fetch(url || "tokenizer.json");
-    if (!r.ok) throw new Error(`tokenizer.json HTTP ${r.status}`);
-    return loadTokenizerData(await r.json());
-  }
-  function toLatin1(s) { const b = new TextEncoder().encode(s); let o = ""; for (const x of b) o += String.fromCharCode(x); return o; }
-  function encode(text) {                      // greedy longest match over bytes
-    const s = toLatin1(text), out = [];
-    let i = 0;
-    while (i < s.length) {
-      let m = null;
-      for (let L = Math.min(tok.maxLen, s.length - i); L > 0; L--) {
-        const sub = s.substr(i, L);
-        if (sub in tok.vocab) { m = sub; break; }
-      }
-      if (m === null) { out.push(tok.unk); i++; continue; }
-      out.push(tok.vocab[m]); i += m.length;
-    }
-    return Int32Array.from(out);
-  }
-  function decode(idArr) {
-    let s = "";
-    for (const id of idArr) {
-      const t = tok.ids[id];
-      if (t === undefined || tok.specials.has(t)) continue;
-      s += t;
-    }
-    const bytes = Uint8Array.from([...s].map(c => c.charCodeAt(0)));
-    return new TextDecoder().decode(bytes);
-  }
-  let CORPUS = buildCorpus();
-  let IDS = encode(CORPUS);
-  let DATASET = "built-in corpus";
-
-  // Training text: FineWeb-Edu (10BT sample), HARDCODED as the only dataset.
-  // The serving Space reads random slices of the parquet shards straight off
-  // the HF CDN with range requests (see server.js /data) — no dependency on
-  // the datasets-server rows API, which 503s routinely. Each device pulls its
-  // own random slice (that's data parallelism — batches were always
-  // per-device anyway). Offline or on failure the built-in corpus stays.
-  const DEFAULT_DS = "HuggingFaceFW/fineweb-edu";
-  async function streamDataset() {              // dataset choice removed on purpose
-    const r = await fetch("data");
-    if (!r.ok) throw new Error(`/data HTTP ${r.status}`);
-    const text = (await r.text()).replace(/[^\x20-\x7e\n]/g, " ");
-    if (text.length < 10000) throw new Error("too little text returned");
-    CORPUS = text.slice(0, 500000);
-    IDS = encode(CORPUS);
-    DATASET = `${DEFAULT_DS} · 10BT sample (parquet via this Space)`;
-    return { name: DATASET, chars: CORPUS.length };
-  }
-  const streamFineWebEdu = () => streamDataset();
-  function datasetName() { return DATASET; }
-
-  // ---- verified matmul: block-scaled INT8 through the units ------------------
-  // CUTLASS ex. 67/81 blockwise scaling: per-row activation scales × per-column
-  // weight scales, one exact LUT/DP4A GEMM, dequant (+ optional fused ReLU) in
-  // the kernel epilogue (ex. 12). Replaces the per-tensor 3-pass: same outlier
-  // robustness at one third of the unit ops.
-  async function vmm(Xf, Wf, m, k, n, ctx, relu) {
-    // ctx.audit re-checks random cells of this LIVE GEMM against the units
-    return V.vgemmBlock(Xf, Wf, { m, k, n, batch: 1, relu: !!relu }, ctx.L, ctx.bgemm, ctx.audit);
-  }
-  // CUTLASS ex. 45 (dual GEMM): sibling GEMMs that share the same LEFT operand
-  // run as ONE batched dispatch, and the shared operand is quantized ONCE
-  // instead of once per sibling. Used for the q/k/v projections — same X
-  // (ln1.y), three weights, identical shapes. Bit-identical to three separate
-  // vmm calls: quantizeRows is deterministic (same input -> same int8+scales),
-  // the tiled copies index exactly like separate batch elements, and block
-  // scales are per-row/per-column PER BATCH ELEMENT, so concatenation changes
-  // no scale and no product. The batched kernel is the same exact-gated bgemm
-  // that training already runs, and the live-shape audit still samples it.
-  async function vmmShared3(Xf, Wa, Wb, Wc, m, k, n, ctx) {
-    const x = V.quantizeRows(Xf, m, k);
-    const xq = new Int8Array(3 * m * k), xs = new Float32Array(3 * m);
-    for (let i = 0; i < 3; i++) { xq.set(x.q, i * m * k); xs.set(x.s, i * m); }
-    const wq = new Int8Array(3 * k * n), ws = new Float32Array(3 * n);
-    [Wa, Wb, Wc].forEach((W, i) => { const w = V.quantizeCols(W, k, n); wq.set(w.q, i * k * n); ws.set(w.s, i * n); });
-    const d = { m, k, n, batch: 3 };
-    let out;
-    if (ctx.bgemm) {
-      out = await ctx.bgemm(xq, wq, xs, ws, d);
-      if (ctx.audit && ctx.audit.due()) {
-        const bad = V.auditTile(xq, wq, xs, ws, d, out, ctx.L, ctx.audit.cells);
-        if (bad) ctx.audit.fail(bad);
-      }
-    } else {
-      out = V.bgemmJS(xq, wq, xs, ws, d, ctx.L);
-    }
-    const MN = m * n;
-    return [out.subarray(0, MN), out.subarray(MN, 2 * MN), out.subarray(2 * MN, 3 * MN)];
-  }
-
-  // ---- layernorm (no affine) -------------------------------------------------
-  function lnFwd(x, rows, C) {
-    const y = new Float32Array(rows * C), sig = new Float32Array(rows);
-    for (let r = 0; r < rows; r++) {
-      let mu = 0; for (let j = 0; j < C; j++) mu += x[r * C + j]; mu /= C;
-      let v = 0; for (let j = 0; j < C; j++) { const d = x[r * C + j] - mu; v += d * d; }
-      const s = Math.sqrt(v / C + 1e-5); sig[r] = s;
-      for (let j = 0; j < C; j++) y[r * C + j] = (x[r * C + j] - mu) / s;
-    }
-    return { y, sig };
-  }
-  function lnBwd(dy, y, sig, rows, C) {
-    const dx = new Float32Array(rows * C);
-    for (let r = 0; r < rows; r++) {
-      let mdy = 0, mdyy = 0;
-      for (let j = 0; j < C; j++) { mdy += dy[r * C + j]; mdyy += dy[r * C + j] * y[r * C + j]; }
-      mdy /= C; mdyy /= C;
-      for (let j = 0; j < C; j++) dx[r * C + j] = (dy[r * C + j] - mdy - y[r * C + j] * mdyy) / sig[r];
-    }
-    return dx;
-  }
-
-  // ---- model -----------------------------------------------------------------
-  // cfg: { c: width, t: seq len, b: batch/device, layers, heads, steps, lr }
-  // engine: the Compute backend object ({bgemm} for the fused WebGPU path) or a
-  // legacy matmulInt8 function (Node tests, inference kit) -> CPU LUT mirror
-  function init(cfg, L, engine, audit) {
-    const c = cfg.c, layers = cfg.layers || 2, heads = cfg.heads || 2, hidden = 2 * c;
-    let seed = 100;
-    const mk = (nEl, scale) => { const w = randn(nEl, mulberry32(seed++)); for (let i = 0; i < nEl; i++) w[i] *= scale; return w; };
-    const params = [], names = [];
-    const add = (name, w) => { params.push(w); names.push(name); return w; };
-    const m = {
-      cfg: { ...cfg, layers, heads, hidden, vocab: vocabSize() },
-      ctx: { L, bgemm: (engine && engine.bgemm) || null,
-             att: (engine && engine.att) || null, fgemm: (engine && engine.fgemm) || null,
-             fgemm2: (engine && engine.fgemm2) || null,
-             mlp: (engine && engine.mlp) || null,
-             audit: audit || null, unitBackward: !!cfg.unitBackward },
-      emb: add("emb", mk(vocabSize() * c, 0.08)),
-      pos: add("pos", mk(cfg.t * c, 0.02)),
-      blocks: [], params, names,
-    };
-    for (let l = 0; l < layers; l++)
-      m.blocks.push({
-        Wq: add(`b${l}.Wq`, mk(c * c, 0.08)), Wk: add(`b${l}.Wk`, mk(c * c, 0.08)),
-        Wv: add(`b${l}.Wv`, mk(c * c, 0.08)), Wo: add(`b${l}.Wo`, mk(c * c, 0.08)),
-        W1: add(`b${l}.W1`, mk(c * hidden, 0.08)), W2: add(`b${l}.W2`, mk(hidden * c, 0.08)),
-      });
-    // weight-tied unembedding: logits use embᵀ (no separate Wu). Halves the
-    // vocab-sized parameters — and with a 16k vocab that's ~half of ALL
-    // parameters, so gradients over the wire shrink ~2× too.
-    m.nParams = params.reduce((a, p) => a + p.length, 0);
-    return m;
-  }
-
-  function sampleBatch(cfg) {
-    const { b, t } = cfg;
-    const X = new Int32Array(b * t), Y = new Int32Array(b * t);
-    for (let i = 0; i < b; i++) {
-      const off = Math.floor(Math.random() * (IDS.length - t - 1));
-      for (let j = 0; j < t; j++) { X[i * t + j] = IDS[off + j]; Y[i * t + j] = IDS[off + j + 1]; }
-    }
-    return { X, Y };
-  }
-
-  // ---- head layout helpers ---------------------------------------------------
-  // q/k/v live as BT×C with head h owning columns [h*hd, (h+1)*hd). The backward
-  // wants every head as its own GEMM problem, so gather once into BH×T×hd and
-  // scatter back at the end — one pass each, instead of slicing per head inside
-  // the loop and paying a GPU dispatch per tiny matmul.
-  function gatherHeads(x, B, T, C, heads, hd) {          // BT×C -> BH×T×hd
-    const out = new Float32Array(B * heads * T * hd);
-    for (let bi = 0; bi < B; bi++)
-      for (let h = 0; h < heads; h++) {
-        const bz = bi * heads + h;
-        for (let ti = 0; ti < T; ti++)
-          for (let j = 0; j < hd; j++) out[(bz * T + ti) * hd + j] = x[(bi * T + ti) * C + h * hd + j];
-      }
-    return out;
-  }
-  function scatterHeadsAcc(dst, src, B, T, C, heads, hd) {   // BH×T×hd -> += BT×C
-    for (let bi = 0; bi < B; bi++)
-      for (let h = 0; h < heads; h++) {
-        const bz = bi * heads + h;
-        for (let ti = 0; ti < T; ti++)
-          for (let j = 0; j < hd; j++) dst[(bi * T + ti) * C + h * hd + j] += src[(bz * T + ti) * hd + j];
-      }
-  }
-  function batchedTranspose(x, batch, rows, cols) {      // per-batch rows×cols -> cols×rows
-    const out = new Float32Array(batch * rows * cols);
-    for (let b = 0; b < batch; b++) {
-      const o = b * rows * cols;
-      for (let r = 0; r < rows; r++)
-        for (let c = 0; c < cols; c++) out[o + c * rows + r] = x[o + r * cols + c];
-    }
-    return out;
-  }
-
-  // ---- forward THROUGH the verified units (caches kept for STE backward) -----
-  async function forward(m, X, Y) {
-    const { c: C, t: T, b: B, layers, heads, hidden, vocab } = m.cfg;
-    const BT = B * T, hd = C / heads, ctx = m.ctx;
-    const cache = { X, Y, blocks: [] };
-    let x = new Float32Array(BT * C);
-    for (let i = 0; i < BT; i++) {
-      const id = X[i], tpos = i % T;
-      for (let j = 0; j < C; j++) x[i * C + j] = m.emb[id * C + j] + m.pos[tpos * C + j];
-    }
-    for (let l = 0; l < layers; l++) {
-      const bl = m.blocks[l], cb = { xin: x };
-      const l1 = lnFwd(x, BT, C); cb.ln1 = l1;
-      // q/k/v share the same left operand — one batched dispatch, one quantize
-      // of ln1.y instead of three (CUTLASS ex. 45; see vmmShared3)
-      const [q, k, v] = await vmmShared3(l1.y, bl.Wq, bl.Wk, bl.Wv, BT, C, C, ctx);
-      cb.q = q; cb.k = k; cb.v = v;
-      const scale = 1 / Math.sqrt(hd);
-      // gather-FUSED attention (CUTLASS ex. 36/52): the kernels read q/k/v in
-      // their natural BT×C layout with head-strided indexing and scatter ctx
-      // straight back — no JS gather copies, no kᵀ transpose. All B×H heads in
-      // one dispatch per stage, every product through the verified units.
-      const BH = B * heads, dAtt = { B, T, heads, hd };
-      // per-(token,head) row quantization: the (BT·heads)×hd view IS the buffer
-      const qq = V.quantizeRows(q, BT * heads, hd), kq = V.quantizeRows(k, BT * heads, hd);
-      const sAll = ctx.att ? await ctx.att.scores(qq.q, kq.q, qq.s, kq.s, dAtt)
-                           : V.attScoresJS(qq.q, kq.q, qq.s, kq.s, dAtt, ctx.L);
-      // live-shape audit: the init gate only ever saw four test shapes
-      if (ctx.att && ctx.audit && ctx.audit.due()) {
-        const bad = V.auditAttScores(qq.q, kq.q, qq.s, kq.s, dAtt, sAll, ctx.L, ctx.audit.cells);
-        if (bad) ctx.audit.fail(bad);
-      }
-      const aAll = new Float32Array(BH * T * T);                   // causal softmax
-      for (let bz = 0; bz < BH; bz++) {
-        const so = bz * T * T;
-        for (let ti = 0; ti < T; ti++) {
-          let mx = -1e30;
-          for (let tj = 0; tj <= ti; tj++) mx = Math.max(mx, sAll[so + ti * T + tj] * scale);
-          let z = 0;
-          for (let tj = 0; tj <= ti; tj++) { const e = Math.exp(sAll[so + ti * T + tj] * scale - mx); aAll[so + ti * T + tj] = e; z += e; }
-          for (let tj = 0; tj <= ti; tj++) aAll[so + ti * T + tj] /= z;
-        }
-      }
-      const aq = V.quantizeRows(aAll, BH * T, T);
-      const vq = V.quantizeHeadCols(v, B, T, heads, hd);
-      const ctxOut = ctx.att ? await ctx.att.ctx(aq.q, vq.q, aq.s, vq.s, dAtt)
-                             : V.attCtxJS(aq.q, vq.q, aq.s, vq.s, dAtt, ctx.L);
-      if (ctx.att && ctx.audit && ctx.audit.due()) {
-        const bad = V.auditAttCtx(aq.q, vq.q, aq.s, vq.s, dAtt, ctxOut, ctx.L, ctx.audit.cells);
-        if (bad) ctx.audit.fail(bad);
-      }
-      cb.aAll = aAll;                                  // backward slices heads from q/k/v/aAll
-      cb.ctxOut = ctxOut;
-      const attnOut = await vmm(ctxOut, bl.Wo, BT, C, C, ctx);
-      const x2 = new Float32Array(BT * C);
-      for (let i = 0; i < x2.length; i++) x2[i] = x[i] + attnOut[i];
-      cb.x2 = x2;
-      const l2 = lnFwd(x2, BT, C); cb.ln2 = l2;
-      // CUTLASS ex. 13 + 23: both MLP GEMMs run back-to-back on the GPU. The
-      // intermediate h1 is quantized ON-DEVICE (exact-gated respec — see
-      // vmlpBlock in verified_core.js) and only its per-row absmax (~1KB)
-      // visits JS between the GEMMs; h1 itself comes back solely because the
-      // STE backward needs it. CPU devices run the bit-identical mirror chain.
-      const { h1, out: mlpOut } = await V.vmlpBlock(l2.y, bl.W1, bl.W2,
-        { m: BT, k: C, h: hidden, n: C }, ctx.L, ctx.mlp, ctx.audit);
-      const mask = new Uint8Array(h1.length);
-      for (let i = 0; i < h1.length; i++) if (h1[i] > 0) mask[i] = 1;
-      cb.h1 = h1; cb.mask = mask;
-      x = new Float32Array(BT * C);
-      for (let i = 0; i < x.length; i++) x[i] = x2[i] + mlpOut[i];
-      cache.blocks.push(cb);
-    }
-    const lf = lnFwd(x, BT, C); cache.lnf = lf; cache.xf = x;
-    const logits = await vmm(lf.y, TC.transpose(m.emb, vocab, C), BT, C, vocab, ctx);  // tied: embᵀ
-    // cross-entropy + dlogits
-    let loss = 0;
-    const dlogits = new Float32Array(BT * vocab);
-    for (let i = 0; i < BT; i++) {
-      let mx = -1e30;
-      for (let j = 0; j < vocab; j++) mx = Math.max(mx, logits[i * vocab + j]);
-      let z = 0;
-      for (let j = 0; j < vocab; j++) z += Math.exp(logits[i * vocab + j] - mx);
-      const lz = Math.log(z) + mx;
-      loss += lz - logits[i * vocab + Y[i]];
-      for (let j = 0; j < vocab; j++)
-        dlogits[i * vocab + j] = (Math.exp(logits[i * vocab + j] - lz) - (j === Y[i] ? 1 : 0)) / BT;
-    }
-    loss /= BT;
-    cache.dlogits = dlogits;
-    return { loss, cache, logits };
-  }
-
-  // ---- STE backward (float), mirrors forward exactly --------------------------
-  // The two vocab-sized matmuls run on the split-K f32 GPU kernel when
-  // available (CUTLASS ex. 06) — same float math, off the JS thread.
-  async function backward(m, cache) {
-    const { c: C, t: T, b: B, layers, heads, hidden, vocab } = m.cfg;
-    const BT = B * T, hd = C / heads, tr = TC.transpose;
-    const g = m.params.map(p => new Float32Array(p.length));
-    const gi = Object.fromEntries(m.names.map((n, i) => [n, i]));
-    // Every matmul here goes through `bmm`. With ctx.unitBackward the STE
-    // gradient is computed BY the verified units (block-scaled int8, exact int32
-    // accumulate) instead of in float. STE is a claim about the math — pretend
-    // the quantizer was the identity — not about the datatype that evaluates it,
-    // so the two are orthogonal and this stays a correct STE.
-    const units = !!m.ctx.unitBackward;
-    const bmm = units
-      ? (A, Bm, mm_, k, n) => vmm(A, Bm, mm_, k, n, m.ctx)
-      : async (A, Bm, mm_, k, n) => TC.matmul(A, Bm, mm_, k, n);
-    // batched: all `batch` problems in ONE dispatch (CUTLASS ex. 05/24). The
-    // per-head backward is 4 GEMMs x B x heads of tiny matrices; issued one at a
-    // time the GPU spends all its time on dispatch overhead rather than math.
-    const bmmB = units
-      ? (A, Bm, rows, k, n, batch) =>
-          V.vgemmBlock(A, Bm, { m: rows, k, n, batch }, m.ctx.L, m.ctx.bgemm, m.ctx.audit)
-      : async (A, Bm, rows, k, n, batch) => {
-          const out = new Float32Array(batch * rows * n);
-          for (let bz = 0; bz < batch; bz++)
-            out.set(TC.matmul(A.subarray(bz * rows * k, (bz + 1) * rows * k),
-                              Bm.subarray(bz * k * n, (bz + 1) * k * n), rows, k, n), bz * rows * n);
-          return out;
-        };
-    // tied unembed: logits = lnf @ embᵀ, so the unembedding gradient flows
-    // straight into emb — dlogitsᵀ @ lnf is V×C, emb's own shape
-    let dlnfIn;
-    if (m.ctx.fgemm2 && !units) {
-      // Both GEMMs consume dlogits (BT x vocab, ~17 MB at the 16512 vocab).
-      // fgemm2 uploads it ONCE and runs both on one submit — profiling had
-      // this pair at 55% of the step, over half of it re-uploading the same
-      // operand. Bit-identical to the two separate calls (gated at init).
-      [g[gi.emb], dlnfIn] = await m.ctx.fgemm2(
-        cache.dlogits,
-        cache.lnf.y, { m: vocab, k: BT, n: C, transA: true },
-        m.emb, { m: BT, k: vocab, n: C });                    // split-K shape
-    } else if (m.ctx.fgemm && !units) {
-      [g[gi.emb], dlnfIn] = await Promise.all([
-        m.ctx.fgemm(cache.dlogits, cache.lnf.y, { m: vocab, k: BT, n: C, transA: true }),
-        m.ctx.fgemm(cache.dlogits, m.emb, { m: BT, k: vocab, n: C }),      // split-K shape
-      ]);
-    } else if (units && m.ctx.bgemm) {
-      // units + GPU: the g.emb operand is dlogitsᵀ (vocab×BT, ~4M elements), and
-      // tr() + quantizeRows() is three full passes over it in JS. Quantizing the
-      // COLUMNS of dlogits directly into transposed int8 is one pass and
-      // bit-identical: same |max| scan, same rounds, in the same order — only
-      // the write pattern changes. The GEMM itself still goes through ctx.bgemm
-      // (exact-gated), and the live-shape audit still samples it.
-      const quantizeColsAsRows = (X, rows, cols) => {   // == quantizeRows(tr(X), cols, rows)
-        const q = new Int8Array(cols * rows), s = new Float32Array(cols);
-        for (let c = 0; c < cols; c++) {
-          let mx = 0;
-          for (let r = 0; r < rows; r++) { const a = Math.abs(X[r * cols + c]); if (a > mx) mx = a; }
-          const sc = Math.max(mx / 127, 1e-8); s[c] = sc;
-          for (let r = 0; r < rows; r++) {
-            const v = Math.round(X[r * cols + c] / sc);
-            q[c * rows + r] = v < -128 ? -128 : v > 127 ? 127 : v;
-          }
-        }
-        return { q, s };
-      };
-      const dlq = quantizeColsAsRows(cache.dlogits, BT, vocab);   // dlogitsᵀ quantized, one pass
-      const wq2 = V.quantizeCols(cache.lnf.y, BT, C);
-      const dEmb = { m: vocab, k: BT, n: C, batch: 1 };
-      const [gEmb, dIn] = await Promise.all([
-        m.ctx.bgemm(dlq.q, wq2.q, dlq.s, wq2.s, dEmb),
-        bmm(cache.dlogits, m.emb, BT, vocab, C),
-      ]);
-      if (m.ctx.audit && m.ctx.audit.due()) {
-        const bad = V.auditTile(dlq.q, wq2.q, dlq.s, wq2.s, dEmb, gEmb, m.ctx.L, m.ctx.audit.cells);
-        if (bad) m.ctx.audit.fail(bad);
-      }
-      g[gi.emb] = gEmb; dlnfIn = dIn;
-    } else {
-      // independent GEMMs — overlap them (these are the two vocab-sized calls,
-      // the largest in the whole backward; each is its own round trip)
-      [g[gi.emb], dlnfIn] = await Promise.all([
-        bmm(tr(cache.dlogits, BT, vocab), cache.lnf.y, vocab, BT, C),
-        bmm(cache.dlogits, m.emb, BT, vocab, C),
-      ]);
-    }
-    let dx = lnBwd(dlnfIn, cache.lnf.y, cache.lnf.sig, BT, C);
-    const scale = 1 / Math.sqrt(hd);
-    // concat helper for fusing sibling GEMMs into one batched dispatch
-    const cat = (...arrs) => {
-      const out = new Float32Array(arrs.reduce((a, x) => a + x.length, 0));
-      let o = 0; for (const x of arrs) { out.set(x, o); o += x.length; }
-      return out;
-    };
-    for (let l = layers - 1; l >= 0; l--) {
-      const bl = m.blocks[l], cb = cache.blocks[l];
-      // mlp: x3 = x2 + relu(ln2 @ W1) @ W2
-      // gW2 and dh1 are independent — overlap their dispatches. On GPU each bmm
-      // is a full upload/submit/readback round trip, so sequential awaits leave
-      // the GPU idle between every pair; this is pure latency, not arithmetic,
-      // and each GEMM's int32 accumulation is exact so overlap changes no bit.
-      const dmlpOut = dx;                              // residual passthrough handled below
-      const [gW2, dh1] = await Promise.all([
-        bmm(tr(cb.h1, BT, hidden), dmlpOut, hidden, BT, C),
-        bmm(dmlpOut, tr(bl.W2, hidden, C), BT, C, hidden),
-      ]);
-      g[gi[`b${l}.W2`]] = gW2;
-      for (let i = 0; i < dh1.length; i++) if (!cb.mask[i]) dh1[i] = 0;
-      const [gW1, dln2raw] = await Promise.all([
-        bmm(tr(cb.ln2.y, BT, C), dh1, C, BT, hidden),
-        bmm(dh1, tr(bl.W1, C, hidden), BT, hidden, C),
-      ]);
-      g[gi[`b${l}.W1`]] = gW1;
-      const dln2in = lnBwd(dln2raw, cb.ln2.y, cb.ln2.sig, BT, C);
-      const dx2 = new Float32Array(BT * C);
-      for (let i = 0; i < dx2.length; i++) dx2[i] = dx[i] + dln2in[i];
-      // attention: x2 = xin + (ctxOut @ Wo)
-      const [gWo, dctx] = await Promise.all([
-        bmm(tr(cb.ctxOut, BT, C), dx2, C, BT, C),
-        bmm(dx2, tr(bl.Wo, C, C), BT, C, C),
-      ]);
-      g[gi[`b${l}.Wo`]] = gWo;
-      // gather every head once, then run each stage as ONE batched GEMM over all
-      // B*heads problems: 4 dispatches per layer instead of 4 per head.
-      const BH = B * heads;
-      const qb = gatherHeads(cb.q, B, T, C, heads, hd);
-      const kb = gatherHeads(cb.k, B, T, C, heads, hd);
-      const vb = gatherHeads(cb.v, B, T, C, heads, hd);
-      const dchb = gatherHeads(dctx, B, T, C, heads, hd);
-      const aT = batchedTranspose(cb.aAll, BH, T, T);        // BH×T×T
-      const vT = batchedTranspose(vb, BH, T, hd);            // BH×hd×T
-      const [dvAll, daAll] = await Promise.all([
-        bmmB(aT, dchb, T, T, hd, BH),                        // aᵀ @ dctx
-        bmmB(dchb, vT, T, hd, T, BH),                        // dctx @ vᵀ
-      ]);
-      // softmax backward is elementwise + a causal row reduction: stays in float
-      // (no matrix math here, so nothing for the units to do)
-      const dsAll = new Float32Array(BH * T * T);
-      for (let bz = 0; bz < BH; bz++) {
-        const o = bz * T * T;
-        for (let ti = 0; ti < T; ti++) {
-          let dot = 0;
-          for (let tj = 0; tj <= ti; tj++) dot += daAll[o + ti * T + tj] * cb.aAll[o + ti * T + tj];
-          for (let tj = 0; tj <= ti; tj++)
-            dsAll[o + ti * T + tj] = cb.aAll[o + ti * T + tj] * (daAll[o + ti * T + tj] - dot) * scale;
-        }
-      }
-      const dsT = batchedTranspose(dsAll, BH, T, T);
-      const [dqAll, dkAll] = await Promise.all([
-        bmmB(dsAll, kb, T, T, hd, BH),                       // ds @ k
-        bmmB(dsT, qb, T, T, hd, BH),                         // dsᵀ @ q
-      ]);
-      const dq = new Float32Array(BT * C), dk = new Float32Array(BT * C), dv = new Float32Array(BT * C);
-      scatterHeadsAcc(dq, dqAll, B, T, C, heads, hd);
-      scatterHeadsAcc(dk, dkAll, B, T, C, heads, hd);
-      scatterHeadsAcc(dv, dvAll, B, T, C, heads, hd);
-      // The QKV weight grads share the same left operand (ln1ᵀ), and the three
-      // dln1in terms share one shape — each trio fuses into ONE batched GEMM
-      // (batch=3) instead of three dispatches. Bit-identical to separate calls:
-      // block scales are per-row of X and per-column of W PER BATCH ELEMENT, so
-      // concatenation changes no scale and no product.
-      const ln1T = tr(cb.ln1.y, BT, C);
-      const [gQKV, dIn3] = await Promise.all([
-        bmmB(cat(ln1T, ln1T, ln1T), cat(dq, dk, dv), C, BT, C, 3),
-        bmmB(cat(dq, dk, dv), cat(tr(bl.Wq, C, C), tr(bl.Wk, C, C), tr(bl.Wv, C, C)), BT, C, C, 3),
-      ]);
-      const CC = C * C, BTC = BT * C;
-      g[gi[`b${l}.Wq`]] = gQKV.slice(0, CC);
-      g[gi[`b${l}.Wk`]] = gQKV.slice(CC, 2 * CC);
-      g[gi[`b${l}.Wv`]] = gQKV.slice(2 * CC, 3 * CC);
-      // sum the three dln1in terms in q,k,v order with an f32 round after EACH
-      // add — the old code accumulated into a Float32Array element three times,
-      // which rounds per step; a bare q+k+v here would run in f64 and round
-      // once, a last-ulp difference that forks replicas. (Exactly the epilogue
-      // mirror lesson: match the rounding schedule, not just the values.)
-      const dln1in = new Float32Array(BTC);
-      for (let i = 0; i < BTC; i++)
-        dln1in[i] = Math.fround(Math.fround(dIn3[i] + dIn3[BTC + i]) + dIn3[2 * BTC + i]);
-      const dxin = lnBwd(dln1in, cb.ln1.y, cb.ln1.sig, BT, C);
-      dx = new Float32Array(BT * C);
-      for (let i = 0; i < dx.length; i++) dx[i] = dx2[i] + dxin[i];
-    }
-    // embedding + positional
-    const ge = g[gi.emb], gp = g[gi.pos];
-    for (let i = 0; i < BT; i++) {
-      const id = cache.X[i], tpos = i % T;
-      for (let j = 0; j < C; j++) { ge[id * C + j] += dx[i * C + j]; gp[tpos * C + j] += dx[i * C + j]; }
-    }
-    // flatten
-    const flat = new Float32Array(m.nParams);
-    let off = 0;
-    for (const t of g) { flat.set(t, off); off += t.length; }
-    return flat;
-  }
-
-  async function trainStep(m) {
-    const { X, Y } = sampleBatch(m.cfg);
-    const { loss, cache } = await forward(m, X, Y);
-    const grad = await backward(m, cache);
-    return { loss, grad };
-  }
-
-  function applyUpdate(m, upd) {          // W -= upd (lr folded in by the optimizer)
-    let off = 0;
-    for (const p of m.params) { for (let i = 0; i < p.length; i++) p[i] -= upd[off + i]; off += p.length; }
-  }
-  function getFlatParams(m) {
-    const flat = new Float32Array(m.nParams);
-    let off = 0;
-    for (const p of m.params) { flat.set(p, off); off += p.length; }
-    return flat;
-  }
-  function setFlatParams(m, flat) {
-    let off = 0;
-    for (const p of m.params) { p.set(flat.subarray(off, off + p.length)); off += p.length; }
-  }
-
-  // greedy sampling — watch the model actually speak
-  async function generate(m, prompt, nChars) {
-    const { t: T } = m.cfg;
-    let ids = [...encode(prompt)];
-    for (let step = 0; step < nChars; step++) {
-      const win = ids.slice(-T);
-      const X = new Int32Array(T), Y = new Int32Array(T);
-      for (let i = 0; i < win.length; i++) X[T - win.length + i] = win[i];
-      const save = m.cfg.b; m.cfg.b = 1;
-      const { logits } = await forward(m, X, Y);
-      m.cfg.b = save;
-      const row = (T - 1) * m.cfg.vocab;
-      let best = 0, bv = -1e30;
-      for (let j = 0; j < m.cfg.vocab; j++) if (logits[row + j] > bv) { bv = logits[row + j]; best = j; }
-      ids.push(best);
-    }
-    return decode(ids);
-  }
-
-  const api = { init, trainStep, applyUpdate, getFlatParams, setFlatParams, generate,
-                streamFineWebEdu, streamDataset, datasetName, loadTokenizer, loadTokenizerData,
-                vocabSize, tokenizerName, encode, decode };
-  if (typeof module !== "undefined" && module.exports) { TC = require("./traincore.js"); V = require("./verified_core.js"); module.exports = api; }
-  else { TC = root.TrainCore; V = root.Verified; root.Transformer = api; }
-})(typeof self !== "undefined" ? self : this);
+// A miniature transformer language model that trains THROUGH the verified INT8
+// units: every matrix product in the forward pass — QKV projections, attention
+// scores, attention·values, output projection, the MLP, and the unembedding —
+// runs through the verified multiply LUT (an emulated INT8 tensor core).
+// Backward is a straight-through estimator in float (the integer path has no
+// gradient), exactly like the Python VerifiedLinear.
+//
+// Task: next-character prediction on a deterministic, self-generated corpus
+// (every peer builds the same text from the same seed — nothing to download).
+(function (root) {
+  "use strict";
+
+  let TC, V;   // TrainCore / Verified — resolved per environment at the end
+
+  function mulberry32(a) { return function () { a |= 0; a = a + 0x6D2B79F5 | 0; let t = Math.imul(a ^ a >>> 15, 1 | a); t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t; return ((t ^ t >>> 14) >>> 0) / 4294967296; }; }
+  function randn(n, rng) { const o = new Float32Array(n); for (let i = 0; i < n; i += 2) { let u = 0, v = 0; while (u === 0) u = rng(); while (v === 0) v = rng(); const m = Math.sqrt(-2 * Math.log(u)); o[i] = m * Math.cos(2 * Math.PI * v); if (i + 1 < n) o[i + 1] = m * Math.sin(2 * Math.PI * v); } return o; }
+
+  // ---- corpus: deterministic cottagecore prose, identical on every peer -----
+  const W_ADJ = ["mossy", "golden", "amber", "quiet", "little", "misty", "sunny", "wild", "cozy", "dusty", "merry", "brave"];
+  const W_NOUN = ["fox", "hare", "owl", "badger", "toad", "sparrow", "otter", "deer", "mushroom", "acorn", "willow", "robin", "river", "meadow", "garden", "lantern"];
+  const W_VERB = ["naps", "sings", "wanders", "hides", "dreams", "waits", "dances", "listens", "rests", "grows"];
+  const W_PREP = ["by", "under", "near", "beside", "beyond", "inside"];
+  function buildCorpus() {
+    const rng = mulberry32(20260712);
+    const pick = (a) => a[Math.floor(rng() * a.length)];
+    let s = "";
+    while (s.length < 60000)
+      s += `the ${pick(W_ADJ)} ${pick(W_NOUN)} ${pick(W_VERB)} ${pick(W_PREP)} the ${pick(W_ADJ)} ${pick(W_NOUN)}. `;
+    return s;
+  }
+  // ---- tokenizer -------------------------------------------------------------
+  // Spikewhale tokenizer (tokenizer.json): byte-level greedy longest-match
+  // ("length-max"), ~16.5k tokens. Until it loads (or if the file is missing)
+  // a 96-char byte-level vocab keeps the app working — but ALL devices in a
+  // group must use the same tokenizer (the config broadcast enforces it).
+  const FALLBACK_CHARS = [...Array(95)].map((_, i) => String.fromCharCode(32 + i)).concat(["\n"]);
+  let tok = {
+    name: "char-96 (fallback)",
+    vocab: Object.fromEntries(FALLBACK_CHARS.map((c, i) => [c, i])),
+    ids: FALLBACK_CHARS, maxLen: 1, size: FALLBACK_CHARS.length,
+    unk: 0, specials: new Set(),
+  };
+  tok.unk = tok.vocab[" "];
+  function vocabSize() { return tok.size; }
+  function tokenizerName() { return tok.name; }
+  function loadTokenizerData(d) {              // plain {vocab, vocab_size, max_token_len}
+    const ids = new Array(d.vocab_size);
+    for (const [t, i] of Object.entries(d.vocab)) ids[i] = t;
+    tok = { name: `Spikewhale length-max (${d.vocab_size} tokens)`,
+            vocab: d.vocab, ids, maxLen: d.max_token_len || 24, size: d.vocab_size,
+            unk: d.vocab["<unk>"] ?? 1,
+            specials: new Set(["<pad>", "<unk>", "<bos>", "<eos>", ...(d.special_tokens || [])]) };
+    IDS = encode(CORPUS);                      // re-tokenize whatever corpus is loaded
+    return tok.name;
+  }
+  async function loadTokenizer(url) {
+    const r = await fetch(url || "tokenizer.json");
+    if (!r.ok) throw new Error(`tokenizer.json HTTP ${r.status}`);
+    return loadTokenizerData(await r.json());
+  }
+  function toLatin1(s) { const b = new TextEncoder().encode(s); let o = ""; for (const x of b) o += String.fromCharCode(x); return o; }
+  function encode(text) {                      // greedy longest match over bytes
+    const s = toLatin1(text), out = [];
+    let i = 0;
+    while (i < s.length) {
+      let m = null;
+      for (let L = Math.min(tok.maxLen, s.length - i); L > 0; L--) {
+        const sub = s.substr(i, L);
+        if (sub in tok.vocab) { m = sub; break; }
+      }
+      if (m === null) { out.push(tok.unk); i++; continue; }
+      out.push(tok.vocab[m]); i += m.length;
+    }
+    return Int32Array.from(out);
+  }
+  function decode(idArr) {
+    let s = "";
+    for (const id of idArr) {
+      const t = tok.ids[id];
+      if (t === undefined || tok.specials.has(t)) continue;
+      s += t;
+    }
+    const bytes = Uint8Array.from([...s].map(c => c.charCodeAt(0)));
+    return new TextDecoder().decode(bytes);
+  }
+  let CORPUS = buildCorpus();
+  let IDS = encode(CORPUS);
+  let DATASET = "built-in corpus";
+
+  // Training text: FineWeb-Edu (10BT sample), HARDCODED as the only dataset.
+  // The serving Space reads random slices of the parquet shards straight off
+  // the HF CDN with range requests (see server.js /data) — no dependency on
+  // the datasets-server rows API, which 503s routinely. Each device pulls its
+  // own random slice (that's data parallelism — batches were always
+  // per-device anyway). Offline or on failure the built-in corpus stays.
+  const DEFAULT_DS = "HuggingFaceFW/fineweb-edu";
+  async function streamDataset() {              // dataset choice removed on purpose
+    const r = await fetch("data");
+    if (!r.ok) throw new Error(`/data HTTP ${r.status}`);
+    const text = (await r.text()).replace(/[^\x20-\x7e\n]/g, " ");
+    if (text.length < 10000) throw new Error("too little text returned");
+    CORPUS = text.slice(0, 500000);
+    IDS = encode(CORPUS);
+    DATASET = `${DEFAULT_DS} · 10BT sample (parquet via this Space)`;
+    return { name: DATASET, chars: CORPUS.length };
+  }
+  const streamFineWebEdu = () => streamDataset();
+  function datasetName() { return DATASET; }
+
+  // ---- verified matmul: block-scaled INT8 through the units ------------------
+  // CUTLASS ex. 67/81 blockwise scaling: per-row activation scales × per-column
+  // weight scales, one exact LUT/DP4A GEMM, dequant (+ optional fused ReLU) in
+  // the kernel epilogue (ex. 12). Replaces the per-tensor 3-pass: same outlier
+  // robustness at one third of the unit ops.
+  async function vmm(Xf, Wf, m, k, n, ctx, relu) {
+    // ctx.audit re-checks random cells of this LIVE GEMM against the units
+    return V.vgemmBlock(Xf, Wf, { m, k, n, batch: 1, relu: !!relu }, ctx.L, ctx.bgemm, ctx.audit);
+  }
+  // CUTLASS ex. 45 (dual GEMM): sibling GEMMs that share the same LEFT operand
+  // run as ONE batched dispatch, and the shared operand is quantized ONCE
+  // instead of once per sibling. Used for the q/k/v projections — same X
+  // (ln1.y), three weights, identical shapes. Bit-identical to three separate
+  // vmm calls: quantizeRows is deterministic (same input -> same int8+scales),
+  // the tiled copies index exactly like separate batch elements, and block
+  // scales are per-row/per-column PER BATCH ELEMENT, so concatenation changes
+  // no scale and no product. The batched kernel is the same exact-gated bgemm
+  // that training already runs, and the live-shape audit still samples it.
+  async function vmmShared3(Xf, Wa, Wb, Wc, m, k, n, ctx) {
+    const x = V.quantizeRows(Xf, m, k);
+    const xq = new Int8Array(3 * m * k), xs = new Float32Array(3 * m);
+    for (let i = 0; i < 3; i++) { xq.set(x.q, i * m * k); xs.set(x.s, i * m); }
+    const wq = new Int8Array(3 * k * n), ws = new Float32Array(3 * n);
+    [Wa, Wb, Wc].forEach((W, i) => { const w = V.quantizeCols(W, k, n); wq.set(w.q, i * k * n); ws.set(w.s, i * n); });
+    const d = { m, k, n, batch: 3 };
+    let out;
+    if (ctx.bgemm) {
+      out = await ctx.bgemm(xq, wq, xs, ws, d);
+      if (ctx.audit && ctx.audit.due()) {
+        const bad = V.auditTile(xq, wq, xs, ws, d, out, ctx.L, ctx.audit.cells);
+        if (bad) ctx.audit.fail(bad);
+      }
+    } else {
+      out = V.bgemmJS(xq, wq, xs, ws, d, ctx.L);
+    }
+    const MN = m * n;
+    return [out.subarray(0, MN), out.subarray(MN, 2 * MN), out.subarray(2 * MN, 3 * MN)];
+  }
+
+  // ---- layernorm (no affine) -------------------------------------------------
+  function lnFwd(x, rows, C) {
+    const y = new Float32Array(rows * C), sig = new Float32Array(rows);
+    for (let r = 0; r < rows; r++) {
+      let mu = 0; for (let j = 0; j < C; j++) mu += x[r * C + j]; mu /= C;
+      let v = 0; for (let j = 0; j < C; j++) { const d = x[r * C + j] - mu; v += d * d; }
+      const s = Math.sqrt(v / C + 1e-5); sig[r] = s;
+      for (let j = 0; j < C; j++) y[r * C + j] = (x[r * C + j] - mu) / s;
+    }
+    return { y, sig };
+  }
+  function lnBwd(dy, y, sig, rows, C) {
+    const dx = new Float32Array(rows * C);
+    for (let r = 0; r < rows; r++) {
+      let mdy = 0, mdyy = 0;
+      for (let j = 0; j < C; j++) { mdy += dy[r * C + j]; mdyy += dy[r * C + j] * y[r * C + j]; }
+      mdy /= C; mdyy /= C;
+      for (let j = 0; j < C; j++) dx[r * C + j] = (dy[r * C + j] - mdy - y[r * C + j] * mdyy) / sig[r];
+    }
+    return dx;
+  }
+
+  // ---- model -----------------------------------------------------------------
+  // cfg: { c: width, t: seq len, b: batch/device, layers, heads, steps, lr }
+  // engine: the Compute backend object ({bgemm} for the fused WebGPU path) or a
+  // legacy matmulInt8 function (Node tests, inference kit) -> CPU LUT mirror
+  function init(cfg, L, engine, audit) {
+    const c = cfg.c, layers = cfg.layers || 2, heads = cfg.heads || 2, hidden = 2 * c;
+    let seed = 100;
+    const mk = (nEl, scale) => { const w = randn(nEl, mulberry32(seed++)); for (let i = 0; i < nEl; i++) w[i] *= scale; return w; };
+    const params = [], names = [];
+    const add = (name, w) => { params.push(w); names.push(name); return w; };
+    const m = {
+      cfg: { ...cfg, layers, heads, hidden, vocab: vocabSize() },
+      ctx: { L, bgemm: (engine && engine.bgemm) || null,
+             att: (engine && engine.att) || null, fgemm: (engine && engine.fgemm) || null,
+             fgemm2: (engine && engine.fgemm2) || null,
+             mlp: (engine && engine.mlp) || null,
+             audit: audit || null, unitBackward: !!cfg.unitBackward },
+      emb: add("emb", mk(vocabSize() * c, 0.08)),
+      pos: add("pos", mk(cfg.t * c, 0.02)),
+      blocks: [], params, names,
+    };
+    for (let l = 0; l < layers; l++)
+      m.blocks.push({
+        Wq: add(`b${l}.Wq`, mk(c * c, 0.08)), Wk: add(`b${l}.Wk`, mk(c * c, 0.08)),
+        Wv: add(`b${l}.Wv`, mk(c * c, 0.08)), Wo: add(`b${l}.Wo`, mk(c * c, 0.08)),
+        W1: add(`b${l}.W1`, mk(c * hidden, 0.08)), W2: add(`b${l}.W2`, mk(hidden * c, 0.08)),
+      });
+    // weight-tied unembedding: logits use embᵀ (no separate Wu). Halves the
+    // vocab-sized parameters — and with a 16k vocab that's ~half of ALL
+    // parameters, so gradients over the wire shrink ~2× too.
+    m.nParams = params.reduce((a, p) => a + p.length, 0);
+    return m;
+  }
+
+  function sampleBatch(cfg) {
+    const { b, t } = cfg;
+    const X = new Int32Array(b * t), Y = new Int32Array(b * t);
+    for (let i = 0; i < b; i++) {
+      const off = Math.floor(Math.random() * (IDS.length - t - 1));
+      for (let j = 0; j < t; j++) { X[i * t + j] = IDS[off + j]; Y[i * t + j] = IDS[off + j + 1]; }
+    }
+    return { X, Y };
+  }
+
+  // ---- head layout helpers ---------------------------------------------------
+  // q/k/v live as BT×C with head h owning columns [h*hd, (h+1)*hd). The backward
+  // wants every head as its own GEMM problem, so gather once into BH×T×hd and
+  // scatter back at the end — one pass each, instead of slicing per head inside
+  // the loop and paying a GPU dispatch per tiny matmul.
+  function gatherHeads(x, B, T, C, heads, hd) {          // BT×C -> BH×T×hd
+    const out = new Float32Array(B * heads * T * hd);
+    for (let bi = 0; bi < B; bi++)
+      for (let h = 0; h < heads; h++) {
+        const bz = bi * heads + h;
+        for (let ti = 0; ti < T; ti++)
+          for (let j = 0; j < hd; j++) out[(bz * T + ti) * hd + j] = x[(bi * T + ti) * C + h * hd + j];
+      }
+    return out;
+  }
+  function scatterHeadsAcc(dst, src, B, T, C, heads, hd) {   // BH×T×hd -> += BT×C
+    for (let bi = 0; bi < B; bi++)
+      for (let h = 0; h < heads; h++) {
+        const bz = bi * heads + h;
+        for (let ti = 0; ti < T; ti++)
+          for (let j = 0; j < hd; j++) dst[(bi * T + ti) * C + h * hd + j] += src[(bz * T + ti) * hd + j];
+      }
+  }
+  function batchedTranspose(x, batch, rows, cols) {      // per-batch rows×cols -> cols×rows
+    const out = new Float32Array(batch * rows * cols);
+    for (let b = 0; b < batch; b++) {
+      const o = b * rows * cols;
+      for (let r = 0; r < rows; r++)
+        for (let c = 0; c < cols; c++) out[o + c * rows + r] = x[o + r * cols + c];
+    }
+    return out;
+  }
+
+  // ---- forward THROUGH the verified units (caches kept for STE backward) -----
+  async function forward(m, X, Y) {
+    const { c: C, t: T, b: B, layers, heads, hidden, vocab } = m.cfg;
+    const BT = B * T, hd = C / heads, ctx = m.ctx;
+    const cache = { X, Y, blocks: [] };
+    let x = new Float32Array(BT * C);
+    for (let i = 0; i < BT; i++) {
+      const id = X[i], tpos = i % T;
+      for (let j = 0; j < C; j++) x[i * C + j] = m.emb[id * C + j] + m.pos[tpos * C + j];
+    }
+    for (let l = 0; l < layers; l++) {
+      const bl = m.blocks[l], cb = { xin: x };
+      const l1 = lnFwd(x, BT, C); cb.ln1 = l1;
+      // q/k/v share the same left operand — one batched dispatch, one quantize
+      // of ln1.y instead of three (CUTLASS ex. 45; see vmmShared3)
+      const [q, k, v] = await vmmShared3(l1.y, bl.Wq, bl.Wk, bl.Wv, BT, C, C, ctx);
+      cb.q = q; cb.k = k; cb.v = v;
+      const scale = 1 / Math.sqrt(hd);
+      // gather-FUSED attention (CUTLASS ex. 36/52): the kernels read q/k/v in
+      // their natural BT×C layout with head-strided indexing and scatter ctx
+      // straight back — no JS gather copies, no kᵀ transpose. All B×H heads in
+      // one dispatch per stage, every product through the verified units.
+      const BH = B * heads, dAtt = { B, T, heads, hd };
+      // per-(token,head) row quantization: the (BT·heads)×hd view IS the buffer
+      const qq = V.quantizeRows(q, BT * heads, hd), kq = V.quantizeRows(k, BT * heads, hd);
+      const sAll = ctx.att ? await ctx.att.scores(qq.q, kq.q, qq.s, kq.s, dAtt)
+                           : V.attScoresJS(qq.q, kq.q, qq.s, kq.s, dAtt, ctx.L);
+      // live-shape audit: the init gate only ever saw four test shapes
+      if (ctx.att && ctx.audit && ctx.audit.due()) {
+        const bad = V.auditAttScores(qq.q, kq.q, qq.s, kq.s, dAtt, sAll, ctx.L, ctx.audit.cells);
+        if (bad) ctx.audit.fail(bad);
+      }
+      const aAll = new Float32Array(BH * T * T);                   // causal softmax
+      for (let bz = 0; bz < BH; bz++) {
+        const so = bz * T * T;
+        for (let ti = 0; ti < T; ti++) {
+          let mx = -1e30;
+          for (let tj = 0; tj <= ti; tj++) mx = Math.max(mx, sAll[so + ti * T + tj] * scale);
+          let z = 0;
+          for (let tj = 0; tj <= ti; tj++) { const e = Math.exp(sAll[so + ti * T + tj] * scale - mx); aAll[so + ti * T + tj] = e; z += e; }
+          for (let tj = 0; tj <= ti; tj++) aAll[so + ti * T + tj] /= z;
+        }
+      }
+      const aq = V.quantizeRows(aAll, BH * T, T);
+      const vq = V.quantizeHeadCols(v, B, T, heads, hd);
+      const ctxOut = ctx.att ? await ctx.att.ctx(aq.q, vq.q, aq.s, vq.s, dAtt)
+                             : V.attCtxJS(aq.q, vq.q, aq.s, vq.s, dAtt, ctx.L);
+      if (ctx.att && ctx.audit && ctx.audit.due()) {
+        const bad = V.auditAttCtx(aq.q, vq.q, aq.s, vq.s, dAtt, ctxOut, ctx.L, ctx.audit.cells);
+        if (bad) ctx.audit.fail(bad);
+      }
+      cb.aAll = aAll;                                  // backward slices heads from q/k/v/aAll
+      cb.ctxOut = ctxOut;
+      const attnOut = await vmm(ctxOut, bl.Wo, BT, C, C, ctx);
+      const x2 = new Float32Array(BT * C);
+      for (let i = 0; i < x2.length; i++) x2[i] = x[i] + attnOut[i];
+      cb.x2 = x2;
+      const l2 = lnFwd(x2, BT, C); cb.ln2 = l2;
+      // CUTLASS ex. 13 + 23: both MLP GEMMs run back-to-back on the GPU. The
+      // intermediate h1 is quantized ON-DEVICE (exact-gated respec — see
+      // vmlpBlock in verified_core.js) and only its per-row absmax (~1KB)
+      // visits JS between the GEMMs; h1 itself comes back solely because the
+      // STE backward needs it. CPU devices run the bit-identical mirror chain.
+      const { h1, out: mlpOut } = await V.vmlpBlock(l2.y, bl.W1, bl.W2,
+        { m: BT, k: C, h: hidden, n: C }, ctx.L, ctx.mlp, ctx.audit);
+      const mask = new Uint8Array(h1.length);
+      for (let i = 0; i < h1.length; i++) if (h1[i] > 0) mask[i] = 1;
+      cb.h1 = h1; cb.mask = mask;
+      x = new Float32Array(BT * C);
+      for (let i = 0; i < x.length; i++) x[i] = x2[i] + mlpOut[i];
+      cache.blocks.push(cb);
+    }
+    const lf = lnFwd(x, BT, C); cache.lnf = lf; cache.xf = x;
+    const logits = await vmm(lf.y, TC.transpose(m.emb, vocab, C), BT, C, vocab, ctx);  // tied: embᵀ
+    // cross-entropy + dlogits
+    let loss = 0;
+    const dlogits = new Float32Array(BT * vocab);
+    for (let i = 0; i < BT; i++) {
+      let mx = -1e30;
+      for (let j = 0; j < vocab; j++) mx = Math.max(mx, logits[i * vocab + j]);
+      let z = 0;
+      for (let j = 0; j < vocab; j++) z += Math.exp(logits[i * vocab + j] - mx);
+      const lz = Math.log(z) + mx;
+      loss += lz - logits[i * vocab + Y[i]];
+      for (let j = 0; j < vocab; j++)
+        dlogits[i * vocab + j] = (Math.exp(logits[i * vocab + j] - lz) - (j === Y[i] ? 1 : 0)) / BT;
+    }
+    loss /= BT;
+    cache.dlogits = dlogits;
+    return { loss, cache, logits };
+  }
+
+  // ---- STE backward (float), mirrors forward exactly --------------------------
+  // The two vocab-sized matmuls run on the split-K f32 GPU kernel when
+  // available (CUTLASS ex. 06) — same float math, off the JS thread.
+  async function backward(m, cache) {
+    const { c: C, t: T, b: B, layers, heads, hidden, vocab } = m.cfg;
+    const BT = B * T, hd = C / heads, tr = TC.transpose;
+    const g = m.params.map(p => new Float32Array(p.length));
+    const gi = Object.fromEntries(m.names.map((n, i) => [n, i]));
+    // Every matmul here goes through `bmm`. With ctx.unitBackward the STE
+    // gradient is computed BY the verified units (block-scaled int8, exact int32
+    // accumulate) instead of in float. STE is a claim about the math — pretend
+    // the quantizer was the identity — not about the datatype that evaluates it,
+    // so the two are orthogonal and this stays a correct STE.
+    const units = !!m.ctx.unitBackward;
+    const bmm = units
+      ? (A, Bm, mm_, k, n) => vmm(A, Bm, mm_, k, n, m.ctx)
+      : async (A, Bm, mm_, k, n) => TC.matmul(A, Bm, mm_, k, n);
+    // batched: all `batch` problems in ONE dispatch (CUTLASS ex. 05/24). The
+    // per-head backward is 4 GEMMs x B x heads of tiny matrices; issued one at a
+    // time the GPU spends all its time on dispatch overhead rather than math.
+    const bmmB = units
+      ? (A, Bm, rows, k, n, batch) =>
+          V.vgemmBlock(A, Bm, { m: rows, k, n, batch }, m.ctx.L, m.ctx.bgemm, m.ctx.audit)
+      : async (A, Bm, rows, k, n, batch) => {
+          const out = new Float32Array(batch * rows * n);
+          for (let bz = 0; bz < batch; bz++)
+            out.set(TC.matmul(A.subarray(bz * rows * k, (bz + 1) * rows * k),
+                              Bm.subarray(bz * k * n, (bz + 1) * k * n), rows, k, n), bz * rows * n);
+          return out;
+        };
+    // tied unembed: logits = lnf @ embᵀ, so the unembedding gradient flows
+    // straight into emb — dlogitsᵀ @ lnf is V×C, emb's own shape
+    let dlnfIn;
+    if (m.ctx.fgemm2 && !units) {
+      // Both GEMMs consume dlogits (BT x vocab, ~17 MB at the 16512 vocab).
+      // fgemm2 uploads it ONCE and runs both on one submit — profiling had
+      // this pair at 55% of the step, over half of it re-uploading the same
+      // operand. Bit-identical to the two separate calls (gated at init).
+      [g[gi.emb], dlnfIn] = await m.ctx.fgemm2(
+        cache.dlogits,
+        cache.lnf.y, { m: vocab, k: BT, n: C, transA: true },
+        m.emb, { m: BT, k: vocab, n: C });                    // split-K shape
+    } else if (m.ctx.fgemm && !units) {
+      [g[gi.emb], dlnfIn] = await Promise.all([
+        m.ctx.fgemm(cache.dlogits, cache.lnf.y, { m: vocab, k: BT, n: C, transA: true }),
+        m.ctx.fgemm(cache.dlogits, m.emb, { m: BT, k: vocab, n: C }),      // split-K shape
+      ]);
+    } else if (units && m.ctx.bgemm) {
+      // units + GPU: the g.emb operand is dlogitsᵀ (vocab×BT, ~4M elements), and
+      // tr() + quantizeRows() is three full passes over it in JS. Quantizing the
+      // COLUMNS of dlogits directly into transposed int8 is one pass and
+      // bit-identical: same |max| scan, same rounds, in the same order — only
+      // the write pattern changes. The GEMM itself still goes through ctx.bgemm
+      // (exact-gated), and the live-shape audit still samples it.
+      const quantizeColsAsRows = (X, rows, cols) => {   // == quantizeRows(tr(X), cols, rows)
+        V.assertFinite(X, "quantizeColsAsRows");
+        const q = new Int8Array(cols * rows), s = new Float32Array(cols);
+        for (let c = 0; c < cols; c++) {
+          let mx = 0;
+          for (let r = 0; r < rows; r++) { const a = Math.abs(X[r * cols + c]); if (a > mx) mx = a; }
+          const sc = Math.max(mx / 127, 1e-8); s[c] = sc;
+          for (let r = 0; r < rows; r++) {
+            const v = Math.round(X[r * cols + c] / sc);
+            q[c * rows + r] = v < -128 ? -128 : v > 127 ? 127 : v;
+          }
+        }
+        return { q, s };
+      };
+      const dlq = quantizeColsAsRows(cache.dlogits, BT, vocab);   // dlogitsᵀ quantized, one pass
+      const wq2 = V.quantizeCols(cache.lnf.y, BT, C);
+      const dEmb = { m: vocab, k: BT, n: C, batch: 1 };
+      const [gEmb, dIn] = await Promise.all([
+        m.ctx.bgemm(dlq.q, wq2.q, dlq.s, wq2.s, dEmb),
+        bmm(cache.dlogits, m.emb, BT, vocab, C),
+      ]);
+      if (m.ctx.audit && m.ctx.audit.due()) {
+        const bad = V.auditTile(dlq.q, wq2.q, dlq.s, wq2.s, dEmb, gEmb, m.ctx.L, m.ctx.audit.cells);
+        if (bad) m.ctx.audit.fail(bad);
+      }
+      g[gi.emb] = gEmb; dlnfIn = dIn;
+    } else {
+      // independent GEMMs — overlap them (these are the two vocab-sized calls,
+      // the largest in the whole backward; each is its own round trip)
+      [g[gi.emb], dlnfIn] = await Promise.all([
+        bmm(tr(cache.dlogits, BT, vocab), cache.lnf.y, vocab, BT, C),
+        bmm(cache.dlogits, m.emb, BT, vocab, C),
+      ]);
+    }
+    let dx = lnBwd(dlnfIn, cache.lnf.y, cache.lnf.sig, BT, C);
+    const scale = 1 / Math.sqrt(hd);
+    // concat helper for fusing sibling GEMMs into one batched dispatch
+    const cat = (...arrs) => {
+      const out = new Float32Array(arrs.reduce((a, x) => a + x.length, 0));
+      let o = 0; for (const x of arrs) { out.set(x, o); o += x.length; }
+      return out;
+    };
+    for (let l = layers - 1; l >= 0; l--) {
+      const bl = m.blocks[l], cb = cache.blocks[l];
+      // mlp: x3 = x2 + relu(ln2 @ W1) @ W2
+      // gW2 and dh1 are independent — overlap their dispatches. On GPU each bmm
+      // is a full upload/submit/readback round trip, so sequential awaits leave
+      // the GPU idle between every pair; this is pure latency, not arithmetic,
+      // and each GEMM's int32 accumulation is exact so overlap changes no bit.
+      const dmlpOut = dx;                              // residual passthrough handled below
+      const [gW2, dh1] = await Promise.all([
+        bmm(tr(cb.h1, BT, hidden), dmlpOut, hidden, BT, C),
+        bmm(dmlpOut, tr(bl.W2, hidden, C), BT, C, hidden),
+      ]);
+      g[gi[`b${l}.W2`]] = gW2;
+      for (let i = 0; i < dh1.length; i++) if (!cb.mask[i]) dh1[i] = 0;
+      const [gW1, dln2raw] = await Promise.all([
+        bmm(tr(cb.ln2.y, BT, C), dh1, C, BT, hidden),
+        bmm(dh1, tr(bl.W1, C, hidden), BT, hidden, C),
+      ]);
+      g[gi[`b${l}.W1`]] = gW1;
+      const dln2in = lnBwd(dln2raw, cb.ln2.y, cb.ln2.sig, BT, C);
+      const dx2 = new Float32Array(BT * C);
+      for (let i = 0; i < dx2.length; i++) dx2[i] = dx[i] + dln2in[i];
+      // attention: x2 = xin + (ctxOut @ Wo)
+      const [gWo, dctx] = await Promise.all([
+        bmm(tr(cb.ctxOut, BT, C), dx2, C, BT, C),
+        bmm(dx2, tr(bl.Wo, C, C), BT, C, C),
+      ]);
+      g[gi[`b${l}.Wo`]] = gWo;
+      // gather every head once, then run each stage as ONE batched GEMM over all
+      // B*heads problems: 4 dispatches per layer instead of 4 per head.
+      const BH = B * heads;
+      const qb = gatherHeads(cb.q, B, T, C, heads, hd);
+      const kb = gatherHeads(cb.k, B, T, C, heads, hd);
+      const vb = gatherHeads(cb.v, B, T, C, heads, hd);
+      const dchb = gatherHeads(dctx, B, T, C, heads, hd);
+      const aT = batchedTranspose(cb.aAll, BH, T, T);        // BH×T×T
+      const vT = batchedTranspose(vb, BH, T, hd);            // BH×hd×T
+      const [dvAll, daAll] = await Promise.all([
+        bmmB(aT, dchb, T, T, hd, BH),                        // aᵀ @ dctx
+        bmmB(dchb, vT, T, hd, T, BH),                        // dctx @ vᵀ
+      ]);
+      // softmax backward is elementwise + a causal row reduction: stays in float
+      // (no matrix math here, so nothing for the units to do)
+      const dsAll = new Float32Array(BH * T * T);
+      for (let bz = 0; bz < BH; bz++) {
+        const o = bz * T * T;
+        for (let ti = 0; ti < T; ti++) {
+          let dot = 0;
+          for (let tj = 0; tj <= ti; tj++) dot += daAll[o + ti * T + tj] * cb.aAll[o + ti * T + tj];
+          for (let tj = 0; tj <= ti; tj++)
+            dsAll[o + ti * T + tj] = cb.aAll[o + ti * T + tj] * (daAll[o + ti * T + tj] - dot) * scale;
+        }
+      }
+      const dsT = batchedTranspose(dsAll, BH, T, T);
+      const [dqAll, dkAll] = await Promise.all([
+        bmmB(dsAll, kb, T, T, hd, BH),                       // ds @ k
+        bmmB(dsT, qb, T, T, hd, BH),                         // dsᵀ @ q
+      ]);
+      const dq = new Float32Array(BT * C), dk = new Float32Array(BT * C), dv = new Float32Array(BT * C);
+      scatterHeadsAcc(dq, dqAll, B, T, C, heads, hd);
+      scatterHeadsAcc(dk, dkAll, B, T, C, heads, hd);
+      scatterHeadsAcc(dv, dvAll, B, T, C, heads, hd);
+      // The QKV weight grads share the same left operand (ln1ᵀ), and the three
+      // dln1in terms share one shape — each trio fuses into ONE batched GEMM
+      // (batch=3) instead of three dispatches. Bit-identical to separate calls:
+      // block scales are per-row of X and per-column of W PER BATCH ELEMENT, so
+      // concatenation changes no scale and no product.
+      const ln1T = tr(cb.ln1.y, BT, C);
+      const [gQKV, dIn3] = await Promise.all([
+        bmmB(cat(ln1T, ln1T, ln1T), cat(dq, dk, dv), C, BT, C, 3),
+        bmmB(cat(dq, dk, dv), cat(tr(bl.Wq, C, C), tr(bl.Wk, C, C), tr(bl.Wv, C, C)), BT, C, C, 3),
+      ]);
+      const CC = C * C, BTC = BT * C;
+      g[gi[`b${l}.Wq`]] = gQKV.slice(0, CC);
+      g[gi[`b${l}.Wk`]] = gQKV.slice(CC, 2 * CC);
+      g[gi[`b${l}.Wv`]] = gQKV.slice(2 * CC, 3 * CC);
+      // sum the three dln1in terms in q,k,v order with an f32 round after EACH
+      // add — the old code accumulated into a Float32Array element three times,
+      // which rounds per step; a bare q+k+v here would run in f64 and round
+      // once, a last-ulp difference that forks replicas. (Exactly the epilogue
+      // mirror lesson: match the rounding schedule, not just the values.)
+      const dln1in = new Float32Array(BTC);
+      for (let i = 0; i < BTC; i++)
+        dln1in[i] = Math.fround(Math.fround(dIn3[i] + dIn3[BTC + i]) + dIn3[2 * BTC + i]);
+      const dxin = lnBwd(dln1in, cb.ln1.y, cb.ln1.sig, BT, C);
+      dx = new Float32Array(BT * C);
+      for (let i = 0; i < dx.length; i++) dx[i] = dx2[i] + dxin[i];
+    }
+    // embedding + positional
+    const ge = g[gi.emb], gp = g[gi.pos];
+    for (let i = 0; i < BT; i++) {
+      const id = cache.X[i], tpos = i % T;
+      for (let j = 0; j < C; j++) { ge[id * C + j] += dx[i * C + j]; gp[tpos * C + j] += dx[i * C + j]; }
+    }
+    // flatten
+    const flat = new Float32Array(m.nParams);
+    let off = 0;
+    for (const t of g) { flat.set(t, off); off += t.length; }
+    return flat;
+  }
+
+  async function trainStep(m) {
+    const { X, Y } = sampleBatch(m.cfg);
+    const { loss, cache } = await forward(m, X, Y);
+    const grad = await backward(m, cache);
+    return { loss, grad };
+  }
+
+  function applyUpdate(m, upd) {          // W -= upd (lr folded in by the optimizer)
+    let off = 0;
+    for (const p of m.params) { for (let i = 0; i < p.length; i++) p[i] -= upd[off + i]; off += p.length; }
+  }
+  function getFlatParams(m) {
+    const flat = new Float32Array(m.nParams);
+    let off = 0;
+    for (const p of m.params) { flat.set(p, off); off += p.length; }
+    return flat;
+  }
+  function setFlatParams(m, flat) {
+    let off = 0;
+    for (const p of m.params) { p.set(flat.subarray(off, off + p.length)); off += p.length; }
+  }
+
+  // greedy sampling — watch the model actually speak
+  async function generate(m, prompt, nChars) {
+    const { t: T } = m.cfg;
+    let ids = [...encode(prompt)];
+    for (let step = 0; step < nChars; step++) {
+      const win = ids.slice(-T);
+      const X = new Int32Array(T), Y = new Int32Array(T);
+      for (let i = 0; i < win.length; i++) X[T - win.length + i] = win[i];
+      const save = m.cfg.b; m.cfg.b = 1;
+      const { logits } = await forward(m, X, Y);
+      m.cfg.b = save;
+      const row = (T - 1) * m.cfg.vocab;
+      let best = 0, bv = -1e30;
+      for (let j = 0; j < m.cfg.vocab; j++) if (logits[row + j] > bv) { bv = logits[row + j]; best = j; }
+      ids.push(best);
+    }
+    return decode(ids);
+  }
+
+  const api = { init, trainStep, applyUpdate, getFlatParams, setFlatParams, generate,
+                streamFineWebEdu, streamDataset, datasetName, loadTokenizer, loadTokenizerData,
+                vocabSize, tokenizerName, encode, decode };
+  if (typeof module !== "undefined" && module.exports) { TC = require("./traincore.js"); V = require("./verified_core.js"); module.exports = api; }
+  else { TC = root.TrainCore; V = root.Verified; root.Transformer = api; }
+})(typeof self !== "undefined" ? self : this);

web/public/verified_core.js

@@ -1,496 +1,519 @@
-// Verified INT8 compute — the emulated GPU logic, in the browser.
-// A layer's forward runs THROUGH the units: quantize -> LUT multiply -> requant
-// -> optional ReLU -> dequant. Backward is a straight-through estimator (the
-// integer path has no gradient), so ordinary float weights still learn.
-// Same units as the Python/Docker DaisyChain; here they're lookup tables.
-(function (root) {
-  "use strict";
-
-  let TC;   // TrainCore (matmul/transpose) — resolved per environment at the end
-
-  function quantize(X) {
-    let mx = 0; for (let i = 0; i < X.length; i++) { const a = Math.abs(X[i]); if (a > mx) mx = a; }
-    const scale = Math.max(mx / 127, 1e-8);
-    const q = new Int8Array(X.length);
-    for (let i = 0; i < X.length; i++) { let v = Math.round(X[i] / scale); q[i] = v < -128 ? -128 : v > 127 ? 127 : v; }
-    return { q, scale };
-  }
-
-  // int8 matmul via the verified multiply LUT: acc(m×n) = sum_k mulLUT[Xq,Wq]
-  function lutMatmulJS(Xq, Wq, m, k, n, L) {
-    const C = new Int32Array(m * n), mul = L.mul;
-    for (let i = 0; i < m; i++) {
-      for (let p = 0; p < k; p++) {
-        const au = (Xq[i * k + p] & 0xFF) * 256, wo = p * n, co = i * n;
-        for (let j = 0; j < n; j++) C[co + j] += mul[au + (Wq[wo + j] & 0xFF)];
-      }
-    }
-    return C;
-  }
-
-  // ---- 3xINT8 fast-accurate GEMM --------------------------------------------
-  // The CUTLASS example-27 "3xTF32" scheme, ported to the verified units:
-  // split each float into a coarse int8 part plus an int8-quantized residual,
-  // run three EXACT LUT GEMMs (hi·hi, hi·lo, lo·hi), drop the negligible
-  // lo·lo, and recombine. Same big/small decomposition NVIDIA uses to recover
-  // near-fp32 accuracy from TF32 tensor cores — here it recovers ~14-bit
-  // accuracy from the 8-bit units, at 3× the unit ops. Every product still
-  // goes through the verified mul8 LUT.
-  function quantize2(X) {
-    const hi = quantize(X);
-    const r = new Float32Array(X.length);
-    for (let i = 0; i < X.length; i++) r[i] = X[i] - hi.q[i] * hi.scale;
-    const lo = quantize(r);
-    return { hi, lo };
-  }
-  function combine3(hh, hl, lh, x, w, len) {
-    const out = new Float32Array(len);
-    const shh = x.hi.scale * w.hi.scale, shl = x.hi.scale * w.lo.scale, slh = x.lo.scale * w.hi.scale;
-    for (let i = 0; i < len; i++) out[i] = hh[i] * shh + hl[i] * shl + lh[i] * slh;
-    return out;
-  }
-  function lutMatmul3JS(Xf, Wf, m, k, n, L) {          // sync, CPU LUT path
-    const x = quantize2(Xf), w = quantize2(Wf);
-    return combine3(lutMatmulJS(x.hi.q, w.hi.q, m, k, n, L),
-                    lutMatmulJS(x.hi.q, w.lo.q, m, k, n, L),
-                    lutMatmulJS(x.lo.q, w.hi.q, m, k, n, L), x, w, m * n);
-  }
-  async function lutMatmul3(Xf, Wf, m, k, n, L, matmulInt8) {  // any backend
-    const x = quantize2(Xf), w = quantize2(Wf);
-    const mm = matmulInt8 || lutMatmulJS;
-    const [hh, hl, lh] = await Promise.all([
-      mm(x.hi.q, w.hi.q, m, k, n, L),
-      mm(x.hi.q, w.lo.q, m, k, n, L),
-      mm(x.lo.q, w.hi.q, m, k, n, L),
-    ]);
-    return combine3(hh, hl, lh, x, w, m * n);
-  }
-
-  // ---- block-scaled verified GEMM (CUTLASS ex. 67/81 blockwise scaling) ------
-  // Per-ROW scales for the activations and per-COLUMN scales for the weights:
-  // the integer math through the mul8 LUT is completely unchanged — only the
-  // dequant uses rs[row]·cs[col] instead of one tensor-wide product, so a single
-  // outlier no longer crushes the quantization resolution of every other
-  // row/column. One LUT pass at this granularity beats the per-tensor 3-pass.
-  function quantizeRows(X, rows, cols) {
-    const q = new Int8Array(rows * cols), s = new Float32Array(rows);
-    for (let r = 0; r < rows; r++) {
-      let mx = 0;
-      for (let c = 0; c < cols; c++) { const a = Math.abs(X[r * cols + c]); if (a > mx) mx = a; }
-      const sc = Math.max(mx / 127, 1e-8); s[r] = sc;
-      for (let c = 0; c < cols; c++) {
-        const v = Math.round(X[r * cols + c] / sc);
-        q[r * cols + c] = v < -128 ? -128 : v > 127 ? 127 : v;
-      }
-    }
-    return { q, s };
-  }
-  function quantizeCols(W, rows, cols) {
-    const q = new Int8Array(rows * cols), s = new Float32Array(cols);
-    for (let c = 0; c < cols; c++) {
-      let mx = 0;
-      for (let r = 0; r < rows; r++) { const a = Math.abs(W[r * cols + c]); if (a > mx) mx = a; }
-      s[c] = Math.max(mx / 127, 1e-8);
-    }
-    for (let r = 0; r < rows; r++)
-      for (let c = 0; c < cols; c++) {
-        const v = Math.round(W[r * cols + c] / s[c]);
-        q[r * cols + c] = v < -128 ? -128 : v > 127 ? 127 : v;
-      }
-    return { q, s };
-  }
-
-  // ---- B2B MLP chain (CUTLASS ex. 13 two-GEMM fusion + ex. 23 epilogue
-  // reduction), respecced for cross-device exactness ---------------------------
-  // The MLP is the one back-to-back GEMM pair with no layernorm/softmax between
-  // (ReLU is already fused in the epilogue), so the intermediate h1 can be
-  // quantized ON the GPU and fed straight to the second GEMM. Two rules make
-  // that fleet-safe:
-  //   1. The per-row |max| (ex. 23) uses only comparisons — exact on any
-  //      hardware, order-independent — and comes back to JS as ~1KB.
-  //   2. Scale DERIVATION (two divisions) stays in JS f64, which IEEE requires
-  //      to be exactly rounded and is therefore identical on every device.
-  //      WGSL division is only 2.5 ULP — a fork waiting to happen — but WGSL
-  //      multiply/add are correctly rounded and floor/clamp are exact. So the
-  //      quantize step is respecced from  round(x / scale)  to
-  //      floor(f32(x * invScale) + 0.5)  — floor(x+0.5) IS Math.round's tie
-  //      rule — and the fround-stepped mirror below is bit-identical to the
-  //      GPU kernel, which is exact-gated against it at init.
-  // NOTE this changes which int8 a value on a rounding boundary lands on
-  // (≤1 step) vs quantizeRows, so old and new builds cannot co-train — the
-  // divergence guard stops such mixed groups by design.
-  function rowAbsMax(X, rows, cols) {
-    const mx = new Float32Array(rows);
-    for (let r = 0; r < rows; r++) {
-      let m = 0;
-      for (let c = 0; c < cols; c++) { const a = Math.abs(X[r * cols + c]); if (a > m) m = a; }
-      mx[r] = m;
-    }
-    return mx;
-  }
-  function scalesFromAbsMax(mx) {                // f64 divisions: exactly rounded, device-identical
-    const scale = new Float32Array(mx.length), inv = new Float32Array(mx.length);
-    for (let i = 0; i < mx.length; i++) {
-      scale[i] = Math.max(mx[i] / 127, 1e-8);
-      inv[i] = 1 / scale[i];                     // recip of the STORED f32 scale
-    }
-    return { scale, inv };
-  }
-  function quantizeRowsInv(X, rows, cols, inv) { // bit-exact mirror of the GPU quantize kernel
-    const q = new Int8Array(rows * cols);
-    for (let r = 0; r < rows; r++) {
-      const iv = inv[r];
-      for (let c = 0; c < cols; c++) {
-        const n = Math.floor(f32(f32(X[r * cols + c] * iv) + 0.5));
-        q[r * cols + c] = n < -128 ? -128 : n > 127 ? 127 : n;
-      }
-    }
-    return q;
-  }
-  // the chained MLP: X @ W1 -> ReLU (fused) -> absmax -> quantize -> @ W2.
-  //   d = { m, k, h, n };  gpuMlp (from webgpu.js) runs both GEMMs + the
-  // on-GPU quantize with one tiny absmax readback between; without it the CPU
-  // mirror chain runs — SAME math, so mixed GPU/CPU fleets stay bit-identical.
-  async function vmlpBlock(Xf, W1f, W2f, d, L, gpuMlp, audit) {
-    const x = quantizeRows(Xf, d.m, d.k);
-    const w1 = quantizeCols(W1f, d.k, d.h);
-    const w2 = quantizeCols(W2f, d.h, d.n);
-    if (gpuMlp) {
-      const r = await gpuMlp(x.q, w1.q, w2.q, x.s, w1.s, w2.s, d);
-      if (audit && audit.due()) {
-        // audit BOTH live GEMMs: gemm1 against the units directly; gemm2 by
-        // reconstructing its exact operand through the proven quantize mirror
-        const bad1 = auditTile(x.q, w1.q, x.s, w1.s, { m: d.m, k: d.k, n: d.h, relu: true }, r.h1, L, audit.cells);
-        if (bad1) { audit.fail("mlp gemm1: " + bad1); return r; }
-        const sc = scalesFromAbsMax(rowAbsMax(r.h1, d.m, d.h));
-        const hq = quantizeRowsInv(r.h1, d.m, d.h, sc.inv);
-        const bad2 = auditTile(hq, w2.q, sc.scale, w2.s, { m: d.m, k: d.h, n: d.n }, r.out, L, audit.cells);
-        if (bad2) audit.fail("mlp gemm2: " + bad2);
-      }
-      return r;
-    }
-    const h1 = bgemmJS(x.q, w1.q, x.s, w1.s, { m: d.m, k: d.k, n: d.h, batch: 1, relu: true }, L);
-    const sc = scalesFromAbsMax(rowAbsMax(h1, d.m, d.h));
-    const hq = quantizeRowsInv(h1, d.m, d.h, sc.inv);
-    const out = bgemmJS(hq, w2.q, sc.scale, w2.s, { m: d.m, k: d.h, n: d.n, batch: 1 }, L);
-    return { h1, out };
-  }
-
-  // ---- epilogue mirror -------------------------------------------------------
-  // BIT-EXACT mirror of the WGSL epilogue `f32(s) * a * b`. WGSL rounds to f32
-  // after the int->float conversion and after EACH multiply; plain JS would do
-  // the whole chain in f64 and round once, which differs in the last ulp. That
-  // last-ulp gap is what used to force a tolerance into the kernel gates —
-  // mirroring the rounding exactly is what lets the gates compare with `!==`.
-  const f32 = Math.fround;
-  function epi(s, a, b) { return f32(f32(f32(s) * a) * b); }
-
-  // f32 equality at the BIT level: `!==` says -0 === 0, but replicas are
-  // compared by hashing raw bytes, so an audit that can't see the sign of zero
-  // could pass a device that later forks the fleet. (Real ISAs have non-IEEE
-  // modes that flush -0 to +0 — e.g. RDNA2 output modifiers / legacy muls.)
-  const _fb = new Float32Array(1), _ub = new Uint32Array(_fb.buffer);
-  function bitDiff(a, b) { _fb[0] = a; const u = _ub[0]; _fb[0] = b; return u !== _ub[0]; }
-
-  // CPU mirror of the fused GPU kernel: batched int8 GEMM through the LUT with
-  // the epilogue (block dequant + optional ReLU) applied before returning —
-  // exactly what the WGSL kernel does on-device. d.acc=true returns the raw
-  // int32 accumulator instead (the exact oracle the fused kernel normally hides).
-  function bgemmJS(Xq, Wq, rs, cs, d, L) {
-    const { m, k, n } = d, batch = d.batch || 1, relu = !!d.relu, mul = L.mul;
-    const raw = !!d.acc;
-    const out = raw ? new Int32Array(batch * m * n) : new Float32Array(batch * m * n);
-    const acc = new Int32Array(n);
-    for (let bz = 0; bz < batch; bz++) {
-      const xo = bz * m * k, wo = bz * k * n, oo = bz * m * n, co = bz * n;
-      for (let i = 0; i < m; i++) {
-        acc.fill(0);
-        const xrow = xo + i * k;
-        for (let p = 0; p < k; p++) {
-          const au = (Xq[xrow + p] & 0xFF) * 256, wrow = wo + p * n;
-          for (let j = 0; j < n; j++) acc[j] += mul[au + (Wq[wrow + j] & 0xFF)];
-        }
-        const orow = oo + i * n;
-        if (raw) { for (let j = 0; j < n; j++) out[orow + j] = acc[j]; continue; }
-        const rscale = rs[bz * m + i];
-        for (let j = 0; j < n; j++) {
-          const v = epi(acc[j], rscale, cs[co + j]);
-          out[orow + j] = relu && v < 0 ? 0 : v;
-        }
-      }
-    }
-    return out;
-  }
-
-  // Recompute a handful of RANDOM output cells of a live GEMM through the LUT
-  // mirror and compare against what the kernel produced. Sampling cells instead
-  // of whole matrices makes this cheap enough to run continuously, at the real
-  // shapes training uses — not once at boot on toy inputs.
-  // STRATIFIED sampling. Uniformly random cells are the wrong instrument for
-  // the bugs that actually occur here: a bounds-guard off-by-one or a pack-tail
-  // padding bug lives on the LAST row/column, and uniform sampling finds that
-  // with probability ~1/n per cell — at the 16512-wide logits GEMM, never. So
-  // the first cells are the structurally dangerous ones (corners, last row,
-  // last column, last batch) chosen deterministically, and the remainder are
-  // random interior cells that catch diffuse bugs. Same principle as poisoning
-  // the buffer pool: construct the dangerous case, don't wait to land on it.
-  function auditTile(Xq, Wq, rs, cs, d, got, L, nCells) {
-    const { m, k, n } = d, batch = d.batch || 1, relu = !!d.relu, mul = L.mul;
-    const N = nCells || 8;
-    const edges = [[0, m - 1, n - 1], [0, 0, n - 1], [0, m - 1, 0], [0, 0, 0],
-                   [batch - 1, m - 1, n - 1], [batch - 1, 0, 0]];
-    for (let t = 0; t < N; t++) {
-      let bz, i, j;
-      if (t < edges.length) { bz = edges[t][0]; i = edges[t][1]; j = edges[t][2]; }
-      else { bz = (Math.random() * batch) | 0; i = (Math.random() * m) | 0; j = (Math.random() * n) | 0; }
-      let acc = 0;
-      const xrow = bz * m * k + i * k, wo = bz * k * n;
-      for (let p = 0; p < k; p++) acc += mul[(Xq[xrow + p] & 0xFF) * 256 + (Wq[wo + p * n + j] & 0xFF)];
-      let v = epi(acc, rs[bz * m + i], cs[bz * n + j]);
-      if (relu && v < 0) v = 0;
-      const idx = (bz * m + i) * n + j;
-      if (bitDiff(got[idx], v))
-        return `GEMM audit failed at [b${bz},${i},${j}] shape ${m}x${k}x${n}: kernel ${Object.is(got[idx], -0) ? "-0" : got[idx]} vs units ${Object.is(v, -0) ? "-0" : v}`;
-    }
-    return null;
-  }
-
-  // ---- exact mirror of the split-K f32 GEMM ----------------------------------
-  // The f32 backward GEMM was the last kernel gated by a TOLERANCE (allclose at
-  // 1e-3) — and this project's own gate mutation test shows allclose waving
-  // through real bugs. The reason was real though: split-K accumulates in a
-  // different ORDER than a naive reference, so bit-equality against the naive
-  // one is impossible. The fix is the same as the epilogue mirror: reproduce
-  // the kernel's order exactly, then compare with `!==`.
-  //   partials: for z in 0..S-1, sum p in [z*ks, min(k,(z+1)*ks)) in order
-  //   reduce:   sum the S partials in ascending z
-  // `fma` selects the rounding schedule for `s + a*b`: WGSL PERMITS a compiler
-  // to contract that into a fused multiply-add (one rounding) instead of two.
-  // Which one the device does is a fact about the device, so the gate tries
-  // both and reports which matches rather than assuming.
-  function fgemmMirror(A, Bm, d, fma) {
-    const { m, k, n } = d, transA = !!d.transA;
-    const S = k > 4096 ? Math.min(16, Math.ceil(k / 2048)) : 1;
-    const ks = Math.ceil(k / S);
-    const out = new Float32Array(m * n);
-    for (let row = 0; row < m; row++)
-      for (let col = 0; col < n; col++) {
-        let acc = 0;                                    // reduce pass, ascending z
-        for (let z = 0; z < S; z++) {
-          const p0 = z * ks, p1 = Math.min(k, p0 + ks);
-          let s = 0;                                    // one partial, in order
-          for (let p = p0; p < p1; p++) {
-            const a = transA ? A[p * m + row] : A[row * k + p];
-            s = fma ? f32(s + a * Bm[p * n + col])      // single rounding
-                    : f32(s + f32(a * Bm[p * n + col]));
-          }
-          acc = f32(acc + s);
-        }
-        out[row * n + col] = acc;
-      }
-    return out;
-  }
-
-  // ---- live audits for the fused attention kernels ---------------------------
-  // The attention kernels had exact INIT gates but nothing at live shapes —
-  // the exact gap the GEMM audit exists to close, left open on the kernels with
-  // the trickiest indexing (head-strided gather, scatter write-back). These
-  // recompute individual output cells from the units, stratified like
-  // auditTile: last/first token pair, last head, last channel first, then
-  // random. Cost is hd (or T) multiply-adds per cell.
-  function auditAttScores(qq, kq, qs, ks, d, got, L, nCells) {
-    const { B, T, heads, hd } = d, C = heads * hd, mul = L.mul, raw = !!d.acc;
-    const N = nCells || 8;
-    const edges = [[B - 1, heads - 1, T - 1, T - 1], [0, 0, 0, 0],
-                   [0, heads - 1, T - 1, 0], [B - 1, 0, 0, T - 1]];
-    for (let t = 0; t < N; t++) {
-      let bi, h, ti, tj;
-      if (t < edges.length) { bi = edges[t][0]; h = edges[t][1]; ti = edges[t][2]; tj = edges[t][3]; }
-      else { bi = (Math.random() * B) | 0; h = (Math.random() * heads) | 0;
-             ti = (Math.random() * T) | 0; tj = (Math.random() * T) | 0; }
-      const bz = bi * heads + h;
-      const qo = (bi * T + ti) * C + h * hd, ko = (bi * T + tj) * C + h * hd;
-      let acc = 0;
-      for (let p = 0; p < hd; p++) acc += mul[(qq[qo + p] & 0xFF) * 256 + (kq[ko + p] & 0xFF)];
-      const v = raw ? acc : epi(acc, qs[(bi * T + ti) * heads + h], ks[(bi * T + tj) * heads + h]);
-      const idx = (bz * T + ti) * T + tj;
-      if (raw ? got[idx] !== v : bitDiff(got[idx], v))
-        return `att.scores audit failed at [b${bi},h${h},${ti},${tj}] B${B}T${T}H${heads}d${hd}: kernel ${got[idx]} vs units ${v}`;
-    }
-    return null;
-  }
-  function auditAttCtx(aq, vq, as, vs, d, got, L, nCells) {
-    const { B, T, heads, hd } = d, C = heads * hd, mul = L.mul, raw = !!d.acc;
-    const N = nCells || 8;
-    const edges = [[B - 1, heads - 1, T - 1, hd - 1], [0, 0, 0, 0],
-                   [0, heads - 1, T - 1, 0], [B - 1, 0, 0, hd - 1]];
-    for (let t = 0; t < N; t++) {
-      let bi, h, ti, j;
-      if (t < edges.length) { bi = edges[t][0]; h = edges[t][1]; ti = edges[t][2]; j = edges[t][3]; }
-      else { bi = (Math.random() * B) | 0; h = (Math.random() * heads) | 0;
-             ti = (Math.random() * T) | 0; j = (Math.random() * hd) | 0; }
-      const bz = bi * heads + h, ao = (bz * T + ti) * T;
-      let acc = 0;
-      for (let tj = 0; tj < T; tj++)
-        acc += mul[(aq[ao + tj] & 0xFF) * 256 + (vq[(bi * T + tj) * C + h * hd + j] & 0xFF)];
-      const v = raw ? acc : epi(acc, as[bz * T + ti], vs[(bi * heads + h) * hd + j]);
-      const idx = (bi * T + ti) * C + h * hd + j;
-      if (raw ? got[idx] !== v : bitDiff(got[idx], v))
-        return `att.ctx audit failed at [b${bi},h${h},${ti},${j}] B${B}T${T}H${heads}d${hd}: kernel ${got[idx]} vs units ${v}`;
-    }
-    return null;
-  }
-
-  // block-scaled verified GEMM, float in → float out.
-  //   d = { m, k, n, batch=1, relu=false }; X is (batch·m)×k, W is batch×(k×n)
-  // gpuBgemm (from webgpu.js) runs the batched kernel with the fused epilogue;
-  // without it the CPU LUT mirror runs. Every product goes through the units.
-  async function vgemmBlock(Xf, Wf, d, L, gpuBgemm, audit) {
-    const { m, k, n } = d, batch = d.batch || 1;
-    const x = quantizeRows(Xf, batch * m, k);
-    let wq, ws;
-    if (batch === 1) {
-      const w = quantizeCols(Wf, k, n); wq = w.q; ws = w.s;
-    } else {
-      wq = new Int8Array(batch * k * n); ws = new Float32Array(batch * n);
-      for (let bz = 0; bz < batch; bz++) {
-        const w = quantizeCols(Wf.subarray(bz * k * n, (bz + 1) * k * n), k, n);
-        wq.set(w.q, bz * k * n); ws.set(w.s, bz * n);
-      }
-    }
-    if (gpuBgemm) {
-      const out = await gpuBgemm(x.q, wq, x.s, ws, d);
-      // continuous re-verification at LIVE shapes: the boot gate only ever saw
-      // toy inputs, so sample a few real cells against the units as we go
-      if (audit && audit.due()) {
-        const bad = auditTile(x.q, wq, x.s, ws, d, out, L, audit.cells);
-        if (bad) audit.fail(bad);
-      }
-      return out;
-    }
-    return bgemmJS(x.q, wq, x.s, ws, d, L);
-  }
-
-  // ---- gather-fused attention through the units (CUTLASS ex. 36/52) ----------
-  // The kernels read q/k/v/ctx directly in their natural BT×C layout with
-  // head-strided indexing — no JS gather copies, no kᵀ transpose, and the
-  // context write scatters straight back into BT×C. Quantization stays
-  // block-scaled: q/k/a per (token,head) row, v per (head,channel) column.
-  // The (BT·heads)×hd row view of q/k IS the contiguous buffer, so
-  // quantizeRows(q, BT·heads, hd) gives per-(token,head) scales for free.
-  function quantizeHeadCols(v, B, T, heads, hd) {   // per (batch,head,channel) column
-    const C = heads * hd;
-    const q = new Int8Array(B * T * C), s = new Float32Array(B * heads * hd);
-    for (let bi = 0; bi < B; bi++)
-      for (let h = 0; h < heads; h++)
-        for (let j = 0; j < hd; j++) {
-          let mx = 0;
-          for (let ti = 0; ti < T; ti++) {
-            const a = Math.abs(v[(bi * T + ti) * C + h * hd + j]);
-            if (a > mx) mx = a;
-          }
-          const sc = Math.max(mx / 127, 1e-8);
-          s[(bi * heads + h) * hd + j] = sc;
-          for (let ti = 0; ti < T; ti++) {
-            const idx = (bi * T + ti) * C + h * hd + j;
-            const w = Math.round(v[idx] / sc);
-            q[idx] = w < -128 ? -128 : w > 127 ? 127 : w;
-          }
-        }
-    return { q, s };
-  }
-  // scores S[bz,ti,tj] = q_row(bi,ti,h) · k_row(bi,tj,h), every product via the LUT
-  // d.acc=true returns the raw int32 accumulator (exact oracle for the kernel gate)
-  function attScoresJS(qq, kq, qs, ks, d, L) {
-    const { B, T, heads, hd } = d, C = heads * hd, mul = L.mul, raw = !!d.acc;
-    const out = raw ? new Int32Array(B * heads * T * T) : new Float32Array(B * heads * T * T);
-    for (let bi = 0; bi < B; bi++) for (let h = 0; h < heads; h++) {
-      const bz = bi * heads + h;
-      for (let ti = 0; ti < T; ti++) {
-        const qo = (bi * T + ti) * C + h * hd, rscale = qs[(bi * T + ti) * heads + h];
-        for (let tj = 0; tj < T; tj++) {
-          const ko = (bi * T + tj) * C + h * hd;
-          let acc = 0;
-          for (let p = 0; p < hd; p++) acc += mul[(qq[qo + p] & 0xFF) * 256 + (kq[ko + p] & 0xFF)];
-          out[(bz * T + ti) * T + tj] = raw ? acc : epi(acc, rscale, ks[(bi * T + tj) * heads + h]);
-        }
-      }
-    }
-    return out;
-  }
-  // ctx[(bi,ti),(h,j)] = Σ_tj a[bz,ti,tj]·v[(bi,tj),(h,j)] — scatter fused into BT×C
-  function attCtxJS(aq, vq, as, vs, d, L) {
-    const { B, T, heads, hd } = d, C = heads * hd, mul = L.mul, raw = !!d.acc;
-    const out = raw ? new Int32Array(B * T * C) : new Float32Array(B * T * C);
-    for (let bi = 0; bi < B; bi++) for (let h = 0; h < heads; h++) {
-      const bz = bi * heads + h;
-      for (let ti = 0; ti < T; ti++) {
-        const ao = (bz * T + ti) * T, rscale = as[bz * T + ti];
-        for (let j = 0; j < hd; j++) {
-          let acc = 0;
-          for (let tj = 0; tj < T; tj++)
-            acc += mul[(aq[ao + tj] & 0xFF) * 256 + (vq[(bi * T + tj) * C + h * hd + j] & 0xFF)];
-          out[(bi * T + ti) * C + h * hd + j] = raw ? acc : epi(acc, rscale, vs[(bi * heads + h) * hd + j]);
-        }
-      }
-    }
-    return out;
-  }
-
-  // one verified layer forward; returns float out (+ cache for STE backward).
-  // Every product goes through the verified INT8 multiply (mul8 LUT) with exact
-  // int32 accumulation — i.e. an emulated INT8 tensor-core GEMM — then dequant.
-  async function linearFwd(X, W, m, k, n, L, useRelu, matmulInt8) {
-    const xq = quantize(X), wq = quantize(W);
-    const acc = await (matmulInt8 || lutMatmulJS)(xq.q, wq.q, m, k, n, L);  // verified multiply
-    const dq = xq.scale * wq.scale;
-    const out = new Float32Array(m * n);
-    const mask = useRelu ? new Uint8Array(m * n) : null;
-    for (let i = 0; i < m * n; i++) {
-      let v = acc[i] * dq;
-      if (useRelu) { if (v > 0) mask[i] = 1; else v = 0; }
-      out[i] = v;
-    }
-    return { out, mask };
-  }
-
-  // 2-layer MLP: X→H (relu) →dout. Forward through verified units, MSE loss.
-  async function forward(X, y, W1, W2, D, L, matmulInt8) {
-    const { n, din, h, dout } = D;
-    const l1 = await linearFwd(X, W1, n, din, h, L, true, matmulInt8);
-    const l2 = await linearFwd(l1.out, W2, n, h, dout, L, false, matmulInt8);
-    const resid = new Float32Array(n * dout); let loss = 0;
-    for (let i = 0; i < resid.length; i++) { const r = l2.out[i] - y[i]; resid[i] = r; loss += r * r; }
-    loss /= resid.length;
-    return { loss, resid, z1: l1.out, mask1: l1.mask };
-  }
-
-  // STE backward (verified matmul treated as float X@W). Returns flat [gW1, gW2].
-  function backward(X, W1, W2, fwd, D) {
-    const { n, din, h, dout } = D;
-    const { resid, z1, mask1 } = fwd;
-    const s = 2 / n;
-    const dout_ = new Float32Array(resid.length);
-    for (let i = 0; i < resid.length; i++) dout_[i] = resid[i] * s;
-    const mm = TC.matmul, tr = TC.transpose;
-    const gW2 = mm(tr(z1, n, h), dout_, h, n, dout);          // z1ᵀ @ dout
-    const dz1 = mm(dout_, tr(W2, h, dout), n, dout, h);        // dout @ W2ᵀ
-    for (let i = 0; i < dz1.length; i++) if (!mask1[i]) dz1[i] = 0;  // relu grad
-    const gW1 = mm(tr(X, n, din), dz1, din, n, h);            // Xᵀ @ dz1
-    const g = new Float32Array(gW1.length + gW2.length);
-    g.set(gW1, 0); g.set(gW2, gW1.length);
-    return g;
-  }
-
-  function splitApply(W1, W2, gAvg, lr) {
-    for (let i = 0; i < W1.length; i++) W1[i] -= lr * gAvg[i];
-    for (let j = 0; j < W2.length; j++) W2[j] -= lr * gAvg[W1.length + j];
-  }
-
-  const api = { quantize, quantize2, quantizeRows, quantizeCols, quantizeHeadCols, lutMatmulJS, lutMatmul3JS, lutMatmul3,
-                bgemmJS, vgemmBlock, auditTile, epi, attScoresJS, attCtxJS, linearFwd, forward, backward, splitApply,
-                rowAbsMax, scalesFromAbsMax, quantizeRowsInv, vmlpBlock, bitDiff,
-                auditAttScores, auditAttCtx, fgemmMirror };
-  if (typeof module !== "undefined" && module.exports) { TC = require("./traincore.js"); module.exports = api; }
-  else { TC = root.TrainCore; root.Verified = api; }
-})(typeof self !== "undefined" ? self : this);
+// Verified INT8 compute — the emulated GPU logic, in the browser.
+// A layer's forward runs THROUGH the units: quantize -> LUT multiply -> requant
+// -> optional ReLU -> dequant. Backward is a straight-through estimator (the
+// integer path has no gradient), so ordinary float weights still learn.
+// Same units as the Python/Docker DaisyChain; here they're lookup tables.
+(function (root) {
+  "use strict";
+
+  let TC;   // TrainCore (matmul/transpose) — resolved per environment at the end
+
+  // A float -> int8 conversion has no answer for NaN/Inf, and the quantizers below give a
+  // silent wrong one: an Inf makes the |max| scale Infinity and EVERY value quantizes to 0
+  // ([0.5, Inf, 3] -> [0, 0, 0]); a NaN is skipped by `a > mx` and then stored as 0 by the
+  // Int8Array. The Python trainer refuses the same input (qat.py); so does this. One pass,
+  // no effect on finite inputs, so builds with and without it still co-train.
+  function assertFinite(X, where) {
+    for (let i = 0; i < X.length; i++) {
+      if (!Number.isFinite(X[i])) {
+        let bad = 0;
+        for (let j = i; j < X.length; j++) if (!Number.isFinite(X[j])) bad++;
+        const e = new Error("non-finite input to " + where + " (" + bad + " of " + X.length +
+                            " values): the int8 quantize would turn it into 0");
+        e.name = "NonFiniteError";
+        throw e;
+      }
+    }
+  }
+
+  function quantize(X) {
+    assertFinite(X, "quantize");
+    let mx = 0; for (let i = 0; i < X.length; i++) { const a = Math.abs(X[i]); if (a > mx) mx = a; }
+    const scale = Math.max(mx / 127, 1e-8);
+    const q = new Int8Array(X.length);
+    for (let i = 0; i < X.length; i++) { let v = Math.round(X[i] / scale); q[i] = v < -128 ? -128 : v > 127 ? 127 : v; }
+    return { q, scale };
+  }
+
+  // int8 matmul via the verified multiply LUT: acc(m×n) = sum_k mulLUT[Xq,Wq]
+  function lutMatmulJS(Xq, Wq, m, k, n, L) {
+    const C = new Int32Array(m * n), mul = L.mul;
+    for (let i = 0; i < m; i++) {
+      for (let p = 0; p < k; p++) {
+        const au = (Xq[i * k + p] & 0xFF) * 256, wo = p * n, co = i * n;
+        for (let j = 0; j < n; j++) C[co + j] += mul[au + (Wq[wo + j] & 0xFF)];
+      }
+    }
+    return C;
+  }
+
+  // ---- 3xINT8 fast-accurate GEMM --------------------------------------------
+  // The CUTLASS example-27 "3xTF32" scheme, ported to the verified units:
+  // split each float into a coarse int8 part plus an int8-quantized residual,
+  // run three EXACT LUT GEMMs (hi·hi, hi·lo, lo·hi), drop the negligible
+  // lo·lo, and recombine. Same big/small decomposition NVIDIA uses to recover
+  // near-fp32 accuracy from TF32 tensor cores — here it recovers ~14-bit
+  // accuracy from the 8-bit units, at 3× the unit ops. Every product still
+  // goes through the verified mul8 LUT.
+  function quantize2(X) {
+    const hi = quantize(X);
+    const r = new Float32Array(X.length);
+    for (let i = 0; i < X.length; i++) r[i] = X[i] - hi.q[i] * hi.scale;
+    const lo = quantize(r);
+    return { hi, lo };
+  }
+  function combine3(hh, hl, lh, x, w, len) {
+    const out = new Float32Array(len);
+    const shh = x.hi.scale * w.hi.scale, shl = x.hi.scale * w.lo.scale, slh = x.lo.scale * w.hi.scale;
+    for (let i = 0; i < len; i++) out[i] = hh[i] * shh + hl[i] * shl + lh[i] * slh;
+    return out;
+  }
+  function lutMatmul3JS(Xf, Wf, m, k, n, L) {          // sync, CPU LUT path
+    const x = quantize2(Xf), w = quantize2(Wf);
+    return combine3(lutMatmulJS(x.hi.q, w.hi.q, m, k, n, L),
+                    lutMatmulJS(x.hi.q, w.lo.q, m, k, n, L),
+                    lutMatmulJS(x.lo.q, w.hi.q, m, k, n, L), x, w, m * n);
+  }
+  async function lutMatmul3(Xf, Wf, m, k, n, L, matmulInt8) {  // any backend
+    const x = quantize2(Xf), w = quantize2(Wf);
+    const mm = matmulInt8 || lutMatmulJS;
+    const [hh, hl, lh] = await Promise.all([
+      mm(x.hi.q, w.hi.q, m, k, n, L),
+      mm(x.hi.q, w.lo.q, m, k, n, L),
+      mm(x.lo.q, w.hi.q, m, k, n, L),
+    ]);
+    return combine3(hh, hl, lh, x, w, m * n);
+  }
+
+  // ---- block-scaled verified GEMM (CUTLASS ex. 67/81 blockwise scaling) ------
+  // Per-ROW scales for the activations and per-COLUMN scales for the weights:
+  // the integer math through the mul8 LUT is completely unchanged — only the
+  // dequant uses rs[row]·cs[col] instead of one tensor-wide product, so a single
+  // outlier no longer crushes the quantization resolution of every other
+  // row/column. One LUT pass at this granularity beats the per-tensor 3-pass.
+  function quantizeRows(X, rows, cols) {
+    assertFinite(X, "quantizeRows");
+    const q = new Int8Array(rows * cols), s = new Float32Array(rows);
+    for (let r = 0; r < rows; r++) {
+      let mx = 0;
+      for (let c = 0; c < cols; c++) { const a = Math.abs(X[r * cols + c]); if (a > mx) mx = a; }
+      const sc = Math.max(mx / 127, 1e-8); s[r] = sc;
+      for (let c = 0; c < cols; c++) {
+        const v = Math.round(X[r * cols + c] / sc);
+        q[r * cols + c] = v < -128 ? -128 : v > 127 ? 127 : v;
+      }
+    }
+    return { q, s };
+  }
+  function quantizeCols(W, rows, cols) {
+    assertFinite(W, "quantizeCols");
+    const q = new Int8Array(rows * cols), s = new Float32Array(cols);
+    for (let c = 0; c < cols; c++) {
+      let mx = 0;
+      for (let r = 0; r < rows; r++) { const a = Math.abs(W[r * cols + c]); if (a > mx) mx = a; }
+      s[c] = Math.max(mx / 127, 1e-8);
+    }
+    for (let r = 0; r < rows; r++)
+      for (let c = 0; c < cols; c++) {
+        const v = Math.round(W[r * cols + c] / s[c]);
+        q[r * cols + c] = v < -128 ? -128 : v > 127 ? 127 : v;
+      }
+    return { q, s };
+  }
+
+  // ---- B2B MLP chain (CUTLASS ex. 13 two-GEMM fusion + ex. 23 epilogue
+  // reduction), respecced for cross-device exactness ---------------------------
+  // The MLP is the one back-to-back GEMM pair with no layernorm/softmax between
+  // (ReLU is already fused in the epilogue), so the intermediate h1 can be
+  // quantized ON the GPU and fed straight to the second GEMM. Two rules make
+  // that fleet-safe:
+  //   1. The per-row |max| (ex. 23) uses only comparisons — exact on any
+  //      hardware, order-independent — and comes back to JS as ~1KB.
+  //   2. Scale DERIVATION (two divisions) stays in JS f64, which IEEE requires
+  //      to be exactly rounded and is therefore identical on every device.
+  //      WGSL division is only 2.5 ULP — a fork waiting to happen — but WGSL
+  //      multiply/add are correctly rounded and floor/clamp are exact. So the
+  //      quantize step is respecced from  round(x / scale)  to
+  //      floor(f32(x * invScale) + 0.5)  — floor(x+0.5) IS Math.round's tie
+  //      rule — and the fround-stepped mirror below is bit-identical to the
+  //      GPU kernel, which is exact-gated against it at init.
+  // NOTE this changes which int8 a value on a rounding boundary lands on
+  // (≤1 step) vs quantizeRows, so old and new builds cannot co-train — the
+  // divergence guard stops such mixed groups by design.
+  function rowAbsMax(X, rows, cols) {
+    assertFinite(X, "rowAbsMax");
+    const mx = new Float32Array(rows);
+    for (let r = 0; r < rows; r++) {
+      let m = 0;
+      for (let c = 0; c < cols; c++) { const a = Math.abs(X[r * cols + c]); if (a > m) m = a; }
+      mx[r] = m;
+    }
+    return mx;
+  }
+  function scalesFromAbsMax(mx) {                // f64 divisions: exactly rounded, device-identical
+    const scale = new Float32Array(mx.length), inv = new Float32Array(mx.length);
+    for (let i = 0; i < mx.length; i++) {
+      scale[i] = Math.max(mx[i] / 127, 1e-8);
+      inv[i] = 1 / scale[i];                     // recip of the STORED f32 scale
+    }
+    return { scale, inv };
+  }
+  function quantizeRowsInv(X, rows, cols, inv) { // bit-exact mirror of the GPU quantize kernel
+    const q = new Int8Array(rows * cols);
+    for (let r = 0; r < rows; r++) {
+      const iv = inv[r];
+      for (let c = 0; c < cols; c++) {
+        const n = Math.floor(f32(f32(X[r * cols + c] * iv) + 0.5));
+        q[r * cols + c] = n < -128 ? -128 : n > 127 ? 127 : n;
+      }
+    }
+    return q;
+  }
+  // the chained MLP: X @ W1 -> ReLU (fused) -> absmax -> quantize -> @ W2.
+  //   d = { m, k, h, n };  gpuMlp (from webgpu.js) runs both GEMMs + the
+  // on-GPU quantize with one tiny absmax readback between; without it the CPU
+  // mirror chain runs — SAME math, so mixed GPU/CPU fleets stay bit-identical.
+  async function vmlpBlock(Xf, W1f, W2f, d, L, gpuMlp, audit) {
+    const x = quantizeRows(Xf, d.m, d.k);
+    const w1 = quantizeCols(W1f, d.k, d.h);
+    const w2 = quantizeCols(W2f, d.h, d.n);
+    if (gpuMlp) {
+      const r = await gpuMlp(x.q, w1.q, w2.q, x.s, w1.s, w2.s, d);
+      if (audit && audit.due()) {
+        // audit BOTH live GEMMs: gemm1 against the units directly; gemm2 by
+        // reconstructing its exact operand through the proven quantize mirror
+        const bad1 = auditTile(x.q, w1.q, x.s, w1.s, { m: d.m, k: d.k, n: d.h, relu: true }, r.h1, L, audit.cells);
+        if (bad1) { audit.fail("mlp gemm1: " + bad1); return r; }
+        const sc = scalesFromAbsMax(rowAbsMax(r.h1, d.m, d.h));
+        const hq = quantizeRowsInv(r.h1, d.m, d.h, sc.inv);
+        const bad2 = auditTile(hq, w2.q, sc.scale, w2.s, { m: d.m, k: d.h, n: d.n }, r.out, L, audit.cells);
+        if (bad2) audit.fail("mlp gemm2: " + bad2);
+      }
+      return r;
+    }
+    const h1 = bgemmJS(x.q, w1.q, x.s, w1.s, { m: d.m, k: d.k, n: d.h, batch: 1, relu: true }, L);
+    const sc = scalesFromAbsMax(rowAbsMax(h1, d.m, d.h));
+    const hq = quantizeRowsInv(h1, d.m, d.h, sc.inv);
+    const out = bgemmJS(hq, w2.q, sc.scale, w2.s, { m: d.m, k: d.h, n: d.n, batch: 1 }, L);
+    return { h1, out };
+  }
+
+  // ---- epilogue mirror -------------------------------------------------------
+  // BIT-EXACT mirror of the WGSL epilogue `f32(s) * a * b`. WGSL rounds to f32
+  // after the int->float conversion and after EACH multiply; plain JS would do
+  // the whole chain in f64 and round once, which differs in the last ulp. That
+  // last-ulp gap is what used to force a tolerance into the kernel gates —
+  // mirroring the rounding exactly is what lets the gates compare with `!==`.
+  const f32 = Math.fround;
+  function epi(s, a, b) { return f32(f32(f32(s) * a) * b); }
+
+  // f32 equality at the BIT level: `!==` says -0 === 0, but replicas are
+  // compared by hashing raw bytes, so an audit that can't see the sign of zero
+  // could pass a device that later forks the fleet. (Real ISAs have non-IEEE
+  // modes that flush -0 to +0 — e.g. RDNA2 output modifiers / legacy muls.)
+  const _fb = new Float32Array(1), _ub = new Uint32Array(_fb.buffer);
+  function bitDiff(a, b) { _fb[0] = a; const u = _ub[0]; _fb[0] = b; return u !== _ub[0]; }
+
+  // CPU mirror of the fused GPU kernel: batched int8 GEMM through the LUT with
+  // the epilogue (block dequant + optional ReLU) applied before returning —
+  // exactly what the WGSL kernel does on-device. d.acc=true returns the raw
+  // int32 accumulator instead (the exact oracle the fused kernel normally hides).
+  function bgemmJS(Xq, Wq, rs, cs, d, L) {
+    const { m, k, n } = d, batch = d.batch || 1, relu = !!d.relu, mul = L.mul;
+    const raw = !!d.acc;
+    const out = raw ? new Int32Array(batch * m * n) : new Float32Array(batch * m * n);
+    const acc = new Int32Array(n);
+    for (let bz = 0; bz < batch; bz++) {
+      const xo = bz * m * k, wo = bz * k * n, oo = bz * m * n, co = bz * n;
+      for (let i = 0; i < m; i++) {
+        acc.fill(0);
+        const xrow = xo + i * k;
+        for (let p = 0; p < k; p++) {
+          const au = (Xq[xrow + p] & 0xFF) * 256, wrow = wo + p * n;
+          for (let j = 0; j < n; j++) acc[j] += mul[au + (Wq[wrow + j] & 0xFF)];
+        }
+        const orow = oo + i * n;
+        if (raw) { for (let j = 0; j < n; j++) out[orow + j] = acc[j]; continue; }
+        const rscale = rs[bz * m + i];
+        for (let j = 0; j < n; j++) {
+          const v = epi(acc[j], rscale, cs[co + j]);
+          out[orow + j] = relu && v < 0 ? 0 : v;
+        }
+      }
+    }
+    return out;
+  }
+
+  // Recompute a handful of RANDOM output cells of a live GEMM through the LUT
+  // mirror and compare against what the kernel produced. Sampling cells instead
+  // of whole matrices makes this cheap enough to run continuously, at the real
+  // shapes training uses — not once at boot on toy inputs.
+  // STRATIFIED sampling. Uniformly random cells are the wrong instrument for
+  // the bugs that actually occur here: a bounds-guard off-by-one or a pack-tail
+  // padding bug lives on the LAST row/column, and uniform sampling finds that
+  // with probability ~1/n per cell — at the 16512-wide logits GEMM, never. So
+  // the first cells are the structurally dangerous ones (corners, last row,
+  // last column, last batch) chosen deterministically, and the remainder are
+  // random interior cells that catch diffuse bugs. Same principle as poisoning
+  // the buffer pool: construct the dangerous case, don't wait to land on it.
+  function auditTile(Xq, Wq, rs, cs, d, got, L, nCells) {
+    const { m, k, n } = d, batch = d.batch || 1, relu = !!d.relu, mul = L.mul;
+    const N = nCells || 8;
+    const edges = [[0, m - 1, n - 1], [0, 0, n - 1], [0, m - 1, 0], [0, 0, 0],
+                   [batch - 1, m - 1, n - 1], [batch - 1, 0, 0]];
+    for (let t = 0; t < N; t++) {
+      let bz, i, j;
+      if (t < edges.length) { bz = edges[t][0]; i = edges[t][1]; j = edges[t][2]; }
+      else { bz = (Math.random() * batch) | 0; i = (Math.random() * m) | 0; j = (Math.random() * n) | 0; }
+      let acc = 0;
+      const xrow = bz * m * k + i * k, wo = bz * k * n;
+      for (let p = 0; p < k; p++) acc += mul[(Xq[xrow + p] & 0xFF) * 256 + (Wq[wo + p * n + j] & 0xFF)];
+      let v = epi(acc, rs[bz * m + i], cs[bz * n + j]);
+      if (relu && v < 0) v = 0;
+      const idx = (bz * m + i) * n + j;
+      if (bitDiff(got[idx], v))
+        return `GEMM audit failed at [b${bz},${i},${j}] shape ${m}x${k}x${n}: kernel ${Object.is(got[idx], -0) ? "-0" : got[idx]} vs units ${Object.is(v, -0) ? "-0" : v}`;
+    }
+    return null;
+  }
+
+  // ---- exact mirror of the split-K f32 GEMM ----------------------------------
+  // The f32 backward GEMM was the last kernel gated by a TOLERANCE (allclose at
+  // 1e-3) — and this project's own gate mutation test shows allclose waving
+  // through real bugs. The reason was real though: split-K accumulates in a
+  // different ORDER than a naive reference, so bit-equality against the naive
+  // one is impossible. The fix is the same as the epilogue mirror: reproduce
+  // the kernel's order exactly, then compare with `!==`.
+  //   partials: for z in 0..S-1, sum p in [z*ks, min(k,(z+1)*ks)) in order
+  //   reduce:   sum the S partials in ascending z
+  // `fma` selects the rounding schedule for `s + a*b`: WGSL PERMITS a compiler
+  // to contract that into a fused multiply-add (one rounding) instead of two.
+  // Which one the device does is a fact about the device, so the gate tries
+  // both and reports which matches rather than assuming.
+  function fgemmMirror(A, Bm, d, fma) {
+    const { m, k, n } = d, transA = !!d.transA;
+    const S = k > 4096 ? Math.min(16, Math.ceil(k / 2048)) : 1;
+    const ks = Math.ceil(k / S);
+    const out = new Float32Array(m * n);
+    for (let row = 0; row < m; row++)
+      for (let col = 0; col < n; col++) {
+        let acc = 0;                                    // reduce pass, ascending z
+        for (let z = 0; z < S; z++) {
+          const p0 = z * ks, p1 = Math.min(k, p0 + ks);
+          let s = 0;                                    // one partial, in order
+          for (let p = p0; p < p1; p++) {
+            const a = transA ? A[p * m + row] : A[row * k + p];
+            s = fma ? f32(s + a * Bm[p * n + col])      // single rounding
+                    : f32(s + f32(a * Bm[p * n + col]));
+          }
+          acc = f32(acc + s);
+        }
+        out[row * n + col] = acc;
+      }
+    return out;
+  }
+
+  // ---- live audits for the fused attention kernels ---------------------------
+  // The attention kernels had exact INIT gates but nothing at live shapes —
+  // the exact gap the GEMM audit exists to close, left open on the kernels with
+  // the trickiest indexing (head-strided gather, scatter write-back). These
+  // recompute individual output cells from the units, stratified like
+  // auditTile: last/first token pair, last head, last channel first, then
+  // random. Cost is hd (or T) multiply-adds per cell.
+  function auditAttScores(qq, kq, qs, ks, d, got, L, nCells) {
+    const { B, T, heads, hd } = d, C = heads * hd, mul = L.mul, raw = !!d.acc;
+    const N = nCells || 8;
+    const edges = [[B - 1, heads - 1, T - 1, T - 1], [0, 0, 0, 0],
+                   [0, heads - 1, T - 1, 0], [B - 1, 0, 0, T - 1]];
+    for (let t = 0; t < N; t++) {
+      let bi, h, ti, tj;
+      if (t < edges.length) { bi = edges[t][0]; h = edges[t][1]; ti = edges[t][2]; tj = edges[t][3]; }
+      else { bi = (Math.random() * B) | 0; h = (Math.random() * heads) | 0;
+             ti = (Math.random() * T) | 0; tj = (Math.random() * T) | 0; }
+      const bz = bi * heads + h;
+      const qo = (bi * T + ti) * C + h * hd, ko = (bi * T + tj) * C + h * hd;
+      let acc = 0;
+      for (let p = 0; p < hd; p++) acc += mul[(qq[qo + p] & 0xFF) * 256 + (kq[ko + p] & 0xFF)];
+      const v = raw ? acc : epi(acc, qs[(bi * T + ti) * heads + h], ks[(bi * T + tj) * heads + h]);
+      const idx = (bz * T + ti) * T + tj;
+      if (raw ? got[idx] !== v : bitDiff(got[idx], v))
+        return `att.scores audit failed at [b${bi},h${h},${ti},${tj}] B${B}T${T}H${heads}d${hd}: kernel ${got[idx]} vs units ${v}`;
+    }
+    return null;
+  }
+  function auditAttCtx(aq, vq, as, vs, d, got, L, nCells) {
+    const { B, T, heads, hd } = d, C = heads * hd, mul = L.mul, raw = !!d.acc;
+    const N = nCells || 8;
+    const edges = [[B - 1, heads - 1, T - 1, hd - 1], [0, 0, 0, 0],
+                   [0, heads - 1, T - 1, 0], [B - 1, 0, 0, hd - 1]];
+    for (let t = 0; t < N; t++) {
+      let bi, h, ti, j;
+      if (t < edges.length) { bi = edges[t][0]; h = edges[t][1]; ti = edges[t][2]; j = edges[t][3]; }
+      else { bi = (Math.random() * B) | 0; h = (Math.random() * heads) | 0;
+             ti = (Math.random() * T) | 0; j = (Math.random() * hd) | 0; }
+      const bz = bi * heads + h, ao = (bz * T + ti) * T;
+      let acc = 0;
+      for (let tj = 0; tj < T; tj++)
+        acc += mul[(aq[ao + tj] & 0xFF) * 256 + (vq[(bi * T + tj) * C + h * hd + j] & 0xFF)];
+      const v = raw ? acc : epi(acc, as[bz * T + ti], vs[(bi * heads + h) * hd + j]);
+      const idx = (bi * T + ti) * C + h * hd + j;
+      if (raw ? got[idx] !== v : bitDiff(got[idx], v))
+        return `att.ctx audit failed at [b${bi},h${h},${ti},${j}] B${B}T${T}H${heads}d${hd}: kernel ${got[idx]} vs units ${v}`;
+    }
+    return null;
+  }
+
+  // block-scaled verified GEMM, float in → float out.
+  //   d = { m, k, n, batch=1, relu=false }; X is (batch·m)×k, W is batch×(k×n)
+  // gpuBgemm (from webgpu.js) runs the batched kernel with the fused epilogue;
+  // without it the CPU LUT mirror runs. Every product goes through the units.
+  async function vgemmBlock(Xf, Wf, d, L, gpuBgemm, audit) {
+    const { m, k, n } = d, batch = d.batch || 1;
+    const x = quantizeRows(Xf, batch * m, k);
+    let wq, ws;
+    if (batch === 1) {
+      const w = quantizeCols(Wf, k, n); wq = w.q; ws = w.s;
+    } else {
+      wq = new Int8Array(batch * k * n); ws = new Float32Array(batch * n);
+      for (let bz = 0; bz < batch; bz++) {
+        const w = quantizeCols(Wf.subarray(bz * k * n, (bz + 1) * k * n), k, n);
+        wq.set(w.q, bz * k * n); ws.set(w.s, bz * n);
+      }
+    }
+    if (gpuBgemm) {
+      const out = await gpuBgemm(x.q, wq, x.s, ws, d);
+      // continuous re-verification at LIVE shapes: the boot gate only ever saw
+      // toy inputs, so sample a few real cells against the units as we go
+      if (audit && audit.due()) {
+        const bad = auditTile(x.q, wq, x.s, ws, d, out, L, audit.cells);
+        if (bad) audit.fail(bad);
+      }
+      return out;
+    }
+    return bgemmJS(x.q, wq, x.s, ws, d, L);
+  }
+
+  // ---- gather-fused attention through the units (CUTLASS ex. 36/52) ----------
+  // The kernels read q/k/v/ctx directly in their natural BT×C layout with
+  // head-strided indexing — no JS gather copies, no kᵀ transpose, and the
+  // context write scatters straight back into BT×C. Quantization stays
+  // block-scaled: q/k/a per (token,head) row, v per (head,channel) column.
+  // The (BT·heads)×hd row view of q/k IS the contiguous buffer, so
+  // quantizeRows(q, BT·heads, hd) gives per-(token,head) scales for free.
+  function quantizeHeadCols(v, B, T, heads, hd) {   // per (batch,head,channel) column
+    assertFinite(v, "quantizeHeadCols");
+    const C = heads * hd;
+    const q = new Int8Array(B * T * C), s = new Float32Array(B * heads * hd);
+    for (let bi = 0; bi < B; bi++)
+      for (let h = 0; h < heads; h++)
+        for (let j = 0; j < hd; j++) {
+          let mx = 0;
+          for (let ti = 0; ti < T; ti++) {
+            const a = Math.abs(v[(bi * T + ti) * C + h * hd + j]);
+            if (a > mx) mx = a;
+          }
+          const sc = Math.max(mx / 127, 1e-8);
+          s[(bi * heads + h) * hd + j] = sc;
+          for (let ti = 0; ti < T; ti++) {
+            const idx = (bi * T + ti) * C + h * hd + j;
+            const w = Math.round(v[idx] / sc);
+            q[idx] = w < -128 ? -128 : w > 127 ? 127 : w;
+          }
+        }
+    return { q, s };
+  }
+  // scores S[bz,ti,tj] = q_row(bi,ti,h) · k_row(bi,tj,h), every product via the LUT
+  // d.acc=true returns the raw int32 accumulator (exact oracle for the kernel gate)
+  function attScoresJS(qq, kq, qs, ks, d, L) {
+    const { B, T, heads, hd } = d, C = heads * hd, mul = L.mul, raw = !!d.acc;
+    const out = raw ? new Int32Array(B * heads * T * T) : new Float32Array(B * heads * T * T);
+    for (let bi = 0; bi < B; bi++) for (let h = 0; h < heads; h++) {
+      const bz = bi * heads + h;
+      for (let ti = 0; ti < T; ti++) {
+        const qo = (bi * T + ti) * C + h * hd, rscale = qs[(bi * T + ti) * heads + h];
+        for (let tj = 0; tj < T; tj++) {
+          const ko = (bi * T + tj) * C + h * hd;
+          let acc = 0;
+          for (let p = 0; p < hd; p++) acc += mul[(qq[qo + p] & 0xFF) * 256 + (kq[ko + p] & 0xFF)];
+          out[(bz * T + ti) * T + tj] = raw ? acc : epi(acc, rscale, ks[(bi * T + tj) * heads + h]);
+        }
+      }
+    }
+    return out;
+  }
+  // ctx[(bi,ti),(h,j)] = Σ_tj a[bz,ti,tj]·v[(bi,tj),(h,j)] — scatter fused into BT×C
+  function attCtxJS(aq, vq, as, vs, d, L) {
+    const { B, T, heads, hd } = d, C = heads * hd, mul = L.mul, raw = !!d.acc;
+    const out = raw ? new Int32Array(B * T * C) : new Float32Array(B * T * C);
+    for (let bi = 0; bi < B; bi++) for (let h = 0; h < heads; h++) {
+      const bz = bi * heads + h;
+      for (let ti = 0; ti < T; ti++) {
+        const ao = (bz * T + ti) * T, rscale = as[bz * T + ti];
+        for (let j = 0; j < hd; j++) {
+          let acc = 0;
+          for (let tj = 0; tj < T; tj++)
+            acc += mul[(aq[ao + tj] & 0xFF) * 256 + (vq[(bi * T + tj) * C + h * hd + j] & 0xFF)];
+          out[(bi * T + ti) * C + h * hd + j] = raw ? acc : epi(acc, rscale, vs[(bi * heads + h) * hd + j]);
+        }
+      }
+    }
+    return out;
+  }
+
+  // one verified layer forward; returns float out (+ cache for STE backward).
+  // Every product goes through the verified INT8 multiply (mul8 LUT) with exact
+  // int32 accumulation — i.e. an emulated INT8 tensor-core GEMM — then dequant.
+  async function linearFwd(X, W, m, k, n, L, useRelu, matmulInt8) {
+    const xq = quantize(X), wq = quantize(W);
+    const acc = await (matmulInt8 || lutMatmulJS)(xq.q, wq.q, m, k, n, L);  // verified multiply
+    const dq = xq.scale * wq.scale;
+    const out = new Float32Array(m * n);
+    const mask = useRelu ? new Uint8Array(m * n) : null;
+    for (let i = 0; i < m * n; i++) {
+      let v = acc[i] * dq;
+      if (useRelu) { if (v > 0) mask[i] = 1; else v = 0; }
+      out[i] = v;
+    }
+    return { out, mask };
+  }
+
+  // 2-layer MLP: X→H (relu) →dout. Forward through verified units, MSE loss.
+  async function forward(X, y, W1, W2, D, L, matmulInt8) {
+    const { n, din, h, dout } = D;
+    const l1 = await linearFwd(X, W1, n, din, h, L, true, matmulInt8);
+    const l2 = await linearFwd(l1.out, W2, n, h, dout, L, false, matmulInt8);
+    const resid = new Float32Array(n * dout); let loss = 0;
+    for (let i = 0; i < resid.length; i++) { const r = l2.out[i] - y[i]; resid[i] = r; loss += r * r; }
+    loss /= resid.length;
+    return { loss, resid, z1: l1.out, mask1: l1.mask };
+  }
+
+  // STE backward (verified matmul treated as float X@W). Returns flat [gW1, gW2].
+  function backward(X, W1, W2, fwd, D) {
+    const { n, din, h, dout } = D;
+    const { resid, z1, mask1 } = fwd;
+    const s = 2 / n;
+    const dout_ = new Float32Array(resid.length);
+    for (let i = 0; i < resid.length; i++) dout_[i] = resid[i] * s;
+    const mm = TC.matmul, tr = TC.transpose;
+    const gW2 = mm(tr(z1, n, h), dout_, h, n, dout);          // z1ᵀ @ dout
+    const dz1 = mm(dout_, tr(W2, h, dout), n, dout, h);        // dout @ W2ᵀ
+    for (let i = 0; i < dz1.length; i++) if (!mask1[i]) dz1[i] = 0;  // relu grad
+    const gW1 = mm(tr(X, n, din), dz1, din, n, h);            // Xᵀ @ dz1
+    const g = new Float32Array(gW1.length + gW2.length);
+    g.set(gW1, 0); g.set(gW2, gW1.length);
+    return g;
+  }
+
+  function splitApply(W1, W2, gAvg, lr) {
+    for (let i = 0; i < W1.length; i++) W1[i] -= lr * gAvg[i];
+    for (let j = 0; j < W2.length; j++) W2[j] -= lr * gAvg[W1.length + j];
+  }
+
+  const api = { assertFinite, quantize, quantize2, quantizeRows, quantizeCols, quantizeHeadCols, lutMatmulJS, lutMatmul3JS, lutMatmul3,
+                bgemmJS, vgemmBlock, auditTile, epi, attScoresJS, attCtxJS, linearFwd, forward, backward, splitApply,
+                rowAbsMax, scalesFromAbsMax, quantizeRowsInv, vmlpBlock, bitDiff,
+                auditAttScores, auditAttCtx, fgemmMirror };
+  if (typeof module !== "undefined" && module.exports) { TC = require("./traincore.js"); module.exports = api; }
+  else { TC = root.TrainCore; root.Verified = api; }
+})(typeof self !== "undefined" ? self : this);

web/test_nonfinite.js

@@ -0,0 +1,25 @@
+// A float -> int8 quantize has no answer for NaN/Inf. Before the guard, an Inf made the |max|
+// scale Infinity and every value quantized to 0 ([0.5, Inf, 3] -> [0, 0, 0]), and a NaN was
+// skipped by the |max| scan and then stored as 0 by the Int8Array -- training continued on
+// zeros with no error. Every quantizer must now refuse, and finite inputs must be unchanged.
+"use strict";
+const V = require("./public/verified_core.js");
+
+let fails = 0;
+function ck(name, cond) { console.log((cond ? "  ok    " : "  FAIL  ") + name); if (!cond) fails++; }
+function throwsNonFinite(fn) { try { fn(); return false; } catch (e) { return e.name === "NonFiniteError"; } }
+
+const bad = [Float32Array.from([0.5, Infinity, 3, -1]), Float32Array.from([0.5, NaN, 3, -1]),
+             Float32Array.from([-Infinity, 1, 2, 3])];
+for (const X of bad) {
+  const tag = Array.from(X).join(",");
+  ck("quantize refuses [" + tag + "]", throwsNonFinite(() => V.quantize(X)));
+  ck("quantizeRows refuses [" + tag + "]", throwsNonFinite(() => V.quantizeRows(X, 2, 2)));
+  ck("quantizeCols refuses [" + tag + "]", throwsNonFinite(() => V.quantizeCols(X, 2, 2)));
+  ck("rowAbsMax refuses [" + tag + "]", throwsNonFinite(() => V.rowAbsMax(X, 2, 2)));
+}
+const q = V.quantize(Float32Array.from([0.5, -1, 3]));
+ck("finite input quantizes as before: [21, -42, 127]", Array.from(q.q).join() === "21,-42,127");
+
+if (fails) { console.log("NONFINITE TEST FAILED: " + fails); process.exit(1); }
+console.log("NONFINITE TEST PASSED — every quantizer refuses NaN/Inf instead of zeroing.");