fix(MAD-223 G.6.g): off-by-one in ml8_fp32_to_e4m3 overflow → +0.40 PPL

ml8.cu::ml8_fp32_to_e4m3 had `e_out >= 15` instead of `> 15` in the
mantissa-round-up overflow path. E4M3 represents e=15, m=0..6 as the
finite values 256, 288, 320, 352, 384, 416, 448 — all valid. Only the
NaN slot (e=15, m=7) is handled separately by the existing m_e4m3==7
guard. The buggy check was prematurely saturating every fp8 activation
in the |x|≈[208..448] post-scale band to ±448, with per-element error
up to 1.75×.

Impact on Cell E (Qwen3.5-4B, llama-perplexity ctx=512 wikitext):
  PPL before: 10.3862  (Δ vs f16 = +0.441)
  PPL after:   9.9849  (Δ vs f16 = +0.040)  ← essentially lossless

ml8-4 now beats Unsloth UD-Q4_K_XL (10.1345, Δ=+0.190) by 0.15 PPL at
comparable size, while running 1.22× faster than f16 on prefill
(842 t/s vs ~692 t/s) and 1.6× smaller (5.4 GB vs 8.5 GB).

Bisect path that surfaced the bug (preserved as diagnostic scripts):
  - check_ml8_gguf_sidecars.py: byte-SHA on centroids/rotation/meta
    against .pt blob — proved sidecars were clean.
  - compare_hip_vs_python_layer0.py: per-stage tensor dump + diff
    (pre-rotation, post-rotation, final). Rotation matched fp32
    epsilon; final diverged → bug downstream of rotation.
  - compare_quant_only.py: feed identical fp32 to HIP's
    ml8_quantize_activations_kernel and torch.float8_e4m3fn cast.
    52 of 40,960 bytes mismatched — all at the e_out=15 boundary.

Env-gated tensor dump hooks (ML8_DUMP=1) added to ml8.cu for future
kernel divergence debugging. Compiled out at runtime when env unset.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

Author: kmbandy

ggml/src/ggml-cuda/ml8.cu

@@ -14,10 +14,56 @@
 
 #include <cstdio>
 #include <cstdlib>
+#include <atomic>
 #include <mutex>
 #include <unordered_map>
 #include <vector>
 
+// G.6.g.C: debug hooks to dump rotation input + ml8_mul_mat output to /tmp
+// for Python-side bit-equivalence comparison. Set env var ML8_DUMP=1 to
+// enable. First-call-only; the static atomics track which dumps have fired.
+namespace {
+std::atomic<bool> g_ml8_dump_rot_done    {false};
+std::atomic<bool> g_ml8_dump_rotdst_done {false};
+std::atomic<bool> g_ml8_dump_mm_done     {false};
+std::atomic<bool> g_ml8_dump_quant_done  {false};
+
+void ml8_dump_u8(const char * path, const uint8_t * d_ptr, size_t n_elems,
+                cudaStream_t stream, int ndim, const int64_t * shape) {
+    std::vector<uint8_t> host(n_elems);
+    cudaMemcpyAsync(host.data(), d_ptr, n_elems, cudaMemcpyDeviceToHost, stream);
+    cudaStreamSynchronize(stream);
+    FILE * f = std::fopen(path, "wb");
+    if (!f) { std::fprintf(stderr, "[ml8-dump] open %s failed\n", path); return; }
+    std::fwrite(&ndim, sizeof(int32_t), 1, f);
+    std::fwrite(shape, sizeof(int64_t), (size_t) ndim, f);
+    std::fwrite(host.data(), 1, n_elems, f);
+    std::fclose(f);
+    std::fprintf(stderr, "[ml8-dump] wrote %s  ndim=%d  n=%zu\n", path, ndim, n_elems);
+}
+
+bool ml8_dump_enabled() {
+    static const bool e = (std::getenv("ML8_DUMP") != nullptr);
+    return e;
+}
+
+void ml8_dump_fp32(const char * path, const float * d_ptr, size_t n_elems,
+                  cudaStream_t stream, int ndim, const int64_t * shape) {
+    std::vector<float> host(n_elems);
+    cudaMemcpyAsync(host.data(), d_ptr, n_elems * sizeof(float),
+                    cudaMemcpyDeviceToHost, stream);
+    cudaStreamSynchronize(stream);
+    FILE * f = std::fopen(path, "wb");
+    if (!f) { std::fprintf(stderr, "[ml8-dump] open %s failed\n", path); return; }
+    // Header: int32 ndim, int64 * shape, then fp32 data
+    std::fwrite(&ndim, sizeof(int32_t), 1, f);
+    std::fwrite(shape, sizeof(int64_t), (size_t) ndim, f);
+    std::fwrite(host.data(), sizeof(float), n_elems, f);
+    std::fclose(f);
+    std::fprintf(stderr, "[ml8-dump] wrote %s  ndim=%d  n=%zu\n", path, ndim, n_elems);
+}
+} // namespace
+
 // On-disk per-block layout: 4-byte fp32 scale, then QK_ML8/2 = 32 packed
 // nibble bytes covering 64 K-elements. sizeof(block_ml8_4) == 36.
 static constexpr int ML8_BLOCK_BYTES   = (int) sizeof(block_ml8_4);
@@ -251,7 +297,12 @@ static __device__ __forceinline__ uint8_t ml8_fp32_to_e4m3(float xv) {
         if (m_e4m3 == 8) {
             m_e4m3 = 0;
             e_out += 1;
-            if (e_out >= 15) {
+            // G.6.g.C BUGFIX (2026-05-26): was `e_out >= 15`, which prematurely
+            // saturated valid e=15, m=0..6 values (256, 288, ..., 448) to ±448.
+            // Only e>15 (= e_real > 8) overflows the E4M3 finite range. The
+            // m=7 NaN slot is handled by the `m_e4m3 == 7` guard below. This
+            // bug cost ~+0.33 PPL on Cell E vs the Python kernel reference.
+            if (e_out > 15) {
                 return (uint8_t)((sign << 7) | (0xFu << 3) | 0x6u);
             }
         }
@@ -438,6 +489,13 @@ void ggml_cuda_op_ml8_mul_mat(
     ggml_cuda_pool_alloc<uint8_t> a_fp8(ctx.pool(),    (size_t) M_pad * (size_t) K);
     ggml_cuda_pool_alloc<float>   a_scale(ctx.pool(), (size_t) M_pad);
 
+    // G.6.g.C: dump pre-quant fp32 activation that the kernel will see.
+    if (ml8_dump_enabled() && !g_ml8_dump_quant_done.load()) {
+        const int64_t shp[2] = { (int64_t) K, (int64_t) M_pad };
+        ml8_dump_fp32("/tmp/ml8_hip_x_prequant.bin", x_src,
+                      (size_t) M_pad * (size_t) K, stream, 2, shp);
+    }
+
     ggml_cuda_ml8_quantize_activations(
         stream,
         x_src,
@@ -446,6 +504,16 @@ void ggml_cuda_op_ml8_mul_mat(
         M_pad,
         K);
 
+    // G.6.g.C: dump fp8 quantized activations + per-row scale on first call.
+    if (ml8_dump_enabled() && !g_ml8_dump_quant_done.exchange(true)) {
+        const int64_t shp_fp8[2]   = { (int64_t) K,     (int64_t) M_pad };
+        const int64_t shp_scale[1] = { (int64_t) M_pad };
+        ml8_dump_u8("/tmp/ml8_hip_a_fp8.bin", a_fp8.get(),
+                    (size_t) M_pad * (size_t) K, stream, 2, shp_fp8);
+        ml8_dump_fp32("/tmp/ml8_hip_a_scale.bin", a_scale.get(),
+                      (size_t) M_pad, stream, 1, shp_scale);
+    }
+
     // ── 4. Allocate bf16 output (M_pad × N) and launch mt_ml8_gemm.
     ggml_cuda_pool_alloc<nv_bfloat16> c_bf16(ctx.pool(), (size_t) M_pad * (size_t) N);
 
@@ -483,6 +551,13 @@ void ggml_cuda_op_ml8_mul_mat(
     GGML_ASSERT(bf16_to_fp32 != nullptr);
     bf16_to_fp32(c_bf16.get(), (float *) dst->data,
                  (size_t) M * (size_t) N, stream);
+
+    // G.6.g.C: dump final mul_mat output on first call.
+    if (ml8_dump_enabled() && !g_ml8_dump_mm_done.exchange(true)) {
+        const int64_t shape[2] = { (int64_t) N, (int64_t) M };
+        ml8_dump_fp32("/tmp/ml8_hip_y_out.bin", (const float *) dst->data,
+                      (size_t) M * (size_t) N, stream, 2, shape);
+    }
 }
 
 // ─────────────────────────────────────────────────────────────────────
@@ -559,6 +634,15 @@ void ggml_cuda_op_ml8_apply_rotation(
     const int n_tokens = (int) x->ne[1];
     const size_t total_elems = (size_t) n_tokens * (size_t) d_dim;
 
+    // G.6.g.C: dump rotation input (pre-rotation activations) on first call.
+    if (ml8_dump_enabled() && !g_ml8_dump_rot_done.exchange(true)) {
+        const int64_t shape[2] = { (int64_t) d_dim, (int64_t) n_tokens };
+        ml8_dump_fp32("/tmp/ml8_hip_x_in.bin", (const float *) x->data,
+                      total_elems, stream, 2, shape);
+    }
+
+    // (rotation kernel runs below; output dump happens after the kernel returns)
+
     // Step 1: copy X into a scratch Z buffer (FWHT is in-place).
     ggml_cuda_pool_alloc<float> z_buf(ctx.pool(), total_elems);
     CUDA_CHECK(cudaMemcpyAsync(z_buf.get(), x->data,
@@ -578,4 +662,11 @@ void ggml_cuda_op_ml8_apply_rotation(
         (float *) dst->data,
         a_dim,
         b_dim);
+
+    // G.6.g.C: dump rotation output (post-FWHT + H_a^T) on first call.
+    if (ml8_dump_enabled() && !g_ml8_dump_rotdst_done.exchange(true)) {
+        const int64_t shape[2] = { (int64_t) d_dim, (int64_t) n_tokens };
+        ml8_dump_fp32("/tmp/ml8_hip_x_rotated.bin", (const float *) dst->data,
+                      total_elems, stream, 2, shape);
+    }
 }

scripts/calibration/check_ml8_gguf_sidecars.py

@@ -0,0 +1,129 @@
+"""Byte-level check: ml8 GGUF sidecars vs source .pt blob.
+
+Reads centroid LUT, rotation H_a, rotation_meta, and ml8 weight bytes
+from a freshly-converted ml8 GGUF and the corresponding cell calibration
+.pt blob. Computes max_abs_diff and SHA256 on each pair.
+
+Run: python3 scripts/calibration/check_ml8_gguf_sidecars.py
+"""
+from __future__ import annotations
+
+import hashlib
+import sys
+from pathlib import Path
+
+import numpy as np
+import torch
+
+sys.path.insert(0, str(Path(__file__).resolve().parent.parent.parent / "gguf-py"))
+import gguf  # noqa: E402
+
+
+GGUF_PATH = Path("/home/kmbandy/models/Qwen3.5-4B-ml8_4-cellE.gguf")
+PT_PATH   = Path("/home/kmbandy/models/cell-e/model_layers_0_mlp_gate_proj.pt")
+BASE      = "blk.0.ffn_gate"
+WEIGHT_NAME = "blk.0.ffn_gate.weight"
+
+
+def sha(b: bytes) -> str:
+    return hashlib.sha256(b).hexdigest()[:16]
+
+
+def main():
+    print(f"GGUF: {GGUF_PATH}")
+    print(f"PT:   {PT_PATH}")
+    print()
+
+    reader = gguf.GGUFReader(str(GGUF_PATH))
+    blob = torch.load(PT_PATH, map_location="cpu", weights_only=False)
+
+    # Map tensor name -> tensor
+    tmap = {t.name: t for t in reader.tensors}
+
+    # ── 1. centroids
+    cent_gguf = tmap[BASE + ".centroids"]
+    cent_gguf_arr = np.asarray(cent_gguf.data)   # uint8 buffer
+    print(f"  GGUF {cent_gguf.name}  shape={tuple(cent_gguf.shape)}  "
+          f"dtype={cent_gguf.tensor_type}  bytes={cent_gguf_arr.nbytes}")
+
+    cent_pt = blob["centroids_per_group"]  # [n_groups_k, 16] fp32
+    cent_pt_fp8 = cent_pt.to(torch.float32).to(torch.float8_e4m3fn)
+    cent_pt_bytes = cent_pt_fp8.contiguous().view(torch.uint8).numpy()
+    print(f"  PT   centroids_per_group → fp8  shape={tuple(cent_pt_fp8.shape)}  "
+          f"bytes={cent_pt_bytes.nbytes}")
+
+    # Compare
+    if cent_gguf_arr.shape == cent_pt_bytes.shape and (cent_gguf_arr == cent_pt_bytes).all():
+        print(f"  ✓ centroids MATCH  sha={sha(cent_gguf_arr.tobytes())}")
+    else:
+        print(f"  ✗ centroids DIFFER")
+        print(f"    GGUF sha = {sha(cent_gguf_arr.tobytes())}")
+        print(f"    PT   sha = {sha(cent_pt_bytes.tobytes())}")
+        # Dequantize both to fp32 and diff
+        g_fp32 = np.asarray(cent_gguf_arr).view(np.uint8)
+        p_fp32 = cent_pt_bytes.flatten()
+        print(f"    GGUF first 16 bytes: {g_fp32[:16].tolist()}")
+        print(f"    PT   first 16 bytes: {p_fp32[:16].tolist()}")
+    print()
+
+    # ── 2. rotation_h_a
+    rot_gguf = tmap.get(BASE + ".rotation_h_a")
+    rot_pt = blob.get("rotation", {}).get("h_a")
+    if rot_gguf is not None and rot_pt is not None:
+        rot_gguf_fp32 = np.asarray(rot_gguf.data).view(np.float32).reshape(tuple(rot_gguf.shape)[::-1])
+        rot_pt_fp32 = rot_pt.to(torch.float32).contiguous().numpy()
+        print(f"  GGUF {rot_gguf.name}  shape={tuple(rot_gguf.shape)}  "
+              f"py-view-shape={rot_gguf_fp32.shape}")
+        print(f"  PT   rotation.h_a  shape={rot_pt_fp32.shape}")
+        # gguf might have reversed dim order
+        if rot_gguf_fp32.shape != rot_pt_fp32.shape:
+            print(f"    shape mismatch! trying transpose...")
+            if rot_gguf_fp32.T.shape == rot_pt_fp32.shape:
+                rot_gguf_fp32 = rot_gguf_fp32.T
+        diff = np.abs(rot_gguf_fp32 - rot_pt_fp32).max() if rot_gguf_fp32.shape == rot_pt_fp32.shape else float("nan")
+        if diff == 0.0:
+            print(f"  ✓ rotation_h_a EXACT MATCH  sha={sha(rot_gguf_fp32.tobytes())}")
+        else:
+            print(f"  ✗ rotation_h_a DIFFERS  max|diff|={diff:.3e}")
+            print(f"    GGUF first 4 values: {rot_gguf_fp32.flatten()[:4]}")
+            print(f"    PT   first 4 values: {rot_pt_fp32.flatten()[:4]}")
+    else:
+        print(f"  (rotation_h_a missing from one side: gguf={rot_gguf is not None}, pt={rot_pt is not None})")
+    print()
+
+    # ── 3. rotation_meta
+    meta_gguf = tmap.get(BASE + ".rotation_meta")
+    if meta_gguf is not None:
+        meta_arr = np.asarray(meta_gguf.data).view(np.int32)
+        print(f"  GGUF {meta_gguf.name}  shape={tuple(meta_gguf.shape)}  values={meta_arr.tolist()}")
+        # Expected: [a_dim, b_dim, in_features, kind_id=1]
+        rot = blob.get("rotation", {})
+        expected = [int(rot.get("a_dim", 0)), int(rot.get("b_dim", 0)),
+                    int(rot.get("in_features", 0)), 1]
+        print(f"  PT   expected:  {expected}")
+        if meta_arr.tolist() == expected:
+            print(f"  ✓ rotation_meta MATCH")
+        else:
+            print(f"  ✗ rotation_meta DIFFERS")
+    print()
+
+    # ── 4. main ml8 weight: block_ml8_4 format
+    # Skip detailed block-by-block for now — if centroids + scales + indices
+    # all match, ml8 weight is right by construction.
+    w_gguf = tmap[WEIGHT_NAME]
+    print(f"  GGUF {w_gguf.name}  shape={tuple(w_gguf.shape)}  "
+          f"dtype={w_gguf.tensor_type}  bytes={np.asarray(w_gguf.data).nbytes}")
+
+    indices = blob["indices"]  # int8 [N, K]
+    scales = blob["scale_per_group"]  # fp32 [N, n_groups_k]
+    print(f"  PT   indices shape={tuple(indices.shape)}  scale shape={tuple(scales.shape)}")
+
+    # Sample first block bytes
+    w_bytes = np.asarray(w_gguf.data)
+    print(f"  GGUF first 8 bytes (block 0): {w_bytes[:8].tolist()}")
+    print(f"  PT   scale[0, 0] = {scales[0, 0].item()}  "
+          f"(if fp32 → bytes: {scales[0, 0].numpy().tobytes().hex()})")
+
+
+if __name__ == "__main__":
+    main()

scripts/calibration/compare_hip_vs_python_layer0.py

@@ -0,0 +1,113 @@
+"""Per-stage bit-equivalence: HIP vs Python ml8 pipeline (G.6.g.C).
+
+Dumps loaded:
+  /tmp/ml8_hip_x_in.bin       — pre-rotation activations (input to ml8_apply_rotation)
+  /tmp/ml8_hip_x_rotated.bin  — post-rotation activations (input to ml8_mul_mat)
+  /tmp/ml8_hip_y_out.bin      — final ml8_mul_mat output (after bf16→fp32)
+
+Tests:
+  Stage A: Python rotation on x_in  ↔  HIP x_rotated     (isolates rotation kernel)
+  Stage B: Python quant+gemm on HIP x_rotated  ↔  HIP y_out  (isolates quant + gemm)
+"""
+from __future__ import annotations
+
+import struct
+import sys
+from pathlib import Path
+
+import numpy as np
+import torch
+
+sys.path.insert(0, str(Path(__file__).resolve().parent))
+from ml8_io import load_ml8_layer as load_blob  # noqa: E402
+from ml8_runtime import ml8_layer_from_blob, ml8_gemm  # noqa: E402
+from kronecker_rotation import KroneckerRotation  # noqa: E402
+
+
+PT_PATH = Path("/home/kmbandy/models/cell-e/model_layers_0_mlp_gate_proj.pt")
+
+
+def read_dump(path: Path) -> np.ndarray:
+    with open(path, "rb") as f:
+        ndim = struct.unpack("<i", f.read(4))[0]
+        shape = list(struct.unpack(f"<{ndim}q", f.read(8 * ndim)))
+        n = 1
+        for s in shape: n *= s
+        data = np.frombuffer(f.read(4 * n), dtype=np.float32).copy()
+    np_shape = tuple(reversed(shape))  # GGML ne[0]=innermost → reverse for numpy
+    return data.reshape(np_shape)
+
+
+def stage_diff(a: np.ndarray, b: np.ndarray, label: str):
+    if a.shape != b.shape:
+        print(f"  {label}: ✗ SHAPE MISMATCH a={a.shape} b={b.shape}")
+        return
+    diff = np.abs(a - b)
+    denom = np.maximum(np.abs(b), 1e-12)
+    rel = diff / denom
+    a_flat = a.reshape(a.shape[0], -1)
+    b_flat = b.reshape(b.shape[0], -1)
+    cos = (a_flat * b_flat).sum(axis=1) / (
+        np.linalg.norm(a_flat, axis=1) * np.linalg.norm(b_flat, axis=1) + 1e-12)
+    print(f"  {label}:")
+    print(f"    max|diff|     = {diff.max():.6e}")
+    print(f"    mean|diff|    = {diff.mean():.6e}")
+    print(f"    max rel diff  = {rel.max():.6e}")
+    print(f"    cos(row) mean = {cos.mean():.6f}   min = {cos.min():.6f}")
+    flat_idx = diff.argmax()
+    r, c = np.unravel_index(flat_idx, diff.shape)
+    print(f"    worst [{r},{c}]: HIP={a[r,c]:.6f}  PY={b[r,c]:.6f}")
+
+
+def main():
+    print("=== HIP dump load ===")
+    x_in_hip      = read_dump(Path("/tmp/ml8_hip_x_in.bin"))
+    x_rotated_hip = read_dump(Path("/tmp/ml8_hip_x_rotated.bin"))
+    y_out_hip     = read_dump(Path("/tmp/ml8_hip_y_out.bin"))
+    M, K = x_in_hip.shape
+    N = y_out_hip.shape[1]
+    print(f"  M={M}  K={K}  N={N}")
+
+    blob = load_blob(PT_PATH)
+    print(f"\n  blob rotation kind: {blob.get('rotation', {}).get('kind')}")
+    device = torch.device("cuda:0")
+
+    # Build the rotation from the blob (same recipe Ml8Linear uses internally)
+    rotation_dict = blob["rotation"]
+    rotation = KroneckerRotation.from_dict(rotation_dict)
+    # rotation.h_a is on CPU by class contract; move to device explicitly.
+    rotation.h_a = rotation.h_a.to(device=device, dtype=torch.float32)
+
+    # ──────────────── Stage A: rotation only ────────────────
+    print("\n=== Stage A: Python rotation ↔ HIP rotation ===")
+    x_in_torch = torch.from_numpy(x_in_hip).contiguous().to(device)
+    x_rotated_py = rotation.forward(x_in_torch).float().cpu().numpy()
+    stage_diff(x_rotated_hip, x_rotated_py, "rotation output (HIP vs Python dense)")
+
+    # ──────────────── Stage B: quant + gemm starting from HIP's rotated ─────────
+    print("\n=== Stage B: Python quant+gemm (input = HIP rotated_x) ↔ HIP final ===")
+    # Pad to multiple of 16
+    pad_m = ((M + 15) // 16) * 16
+    if pad_m != M:
+        pad_rows = pad_m - M
+        x_rot_padded = np.concatenate([x_rotated_hip, np.zeros((pad_rows, K), dtype=np.float32)], axis=0)
+    else:
+        x_rot_padded = x_rotated_hip
+
+    layer = ml8_layer_from_blob(blob, device=device)
+    x_rot_torch = torch.from_numpy(x_rot_padded).contiguous().to(device)
+
+    # Quant: same math as Ml8Linear.forward
+    FP8_MAX = 448.0
+    row_max = x_rot_torch.abs().amax(dim=1, keepdim=True).clamp(min=1e-8)
+    a_scale = (row_max / FP8_MAX).squeeze(1).contiguous()
+    x_fp8 = (x_rot_torch / a_scale.unsqueeze(1)).to(torch.float8_e4m3fn).contiguous()
+
+    with torch.inference_mode():
+        y_bf16 = ml8_gemm(x_fp8, layer, a_scale=a_scale, out_dtype=torch.bfloat16)
+    y_py = y_bf16.float().cpu().numpy()[:M]
+    stage_diff(y_out_hip, y_py, "final output (HIP vs Python quant+gemm)")
+
+
+if __name__ == "__main__":
+    main()