diff --git a/CMakeLists.txt b/CMakeLists.txt
index 190a0f07..8cba44fb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -357,6 +357,16 @@ target_link_libraries(benchmarks libgemma hwy hwy_contrib nlohmann_json::nlohman
 add_executable(debug_prompt evals/debug_prompt.cc)
 target_link_libraries(debug_prompt libgemma hwy hwy_contrib nlohmann_json::nlohmann_json)
 
+add_library(model_comparison evals/model_comparison.cc)
+target_include_directories(model_comparison PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
+
+add_executable(model_comparison_test evals/model_comparison_test.cc)
+target_link_libraries(model_comparison_test model_comparison)
+
+add_executable(gemma_mmlu evals/run_mmlu.cc)
+target_link_libraries(gemma_mmlu libgemma model_comparison hwy hwy_contrib
+                      nlohmann_json::nlohmann_json)
+
 ## Tests
 set(GEMMA_ENABLE_TESTS OFF CACHE BOOL "Enable Gemma tests")
 if (GEMMA_ENABLE_TESTS)
@@ -445,6 +455,27 @@ endif()  # GEMMA_ENABLE_TESTS
 
 ## Tools
 
+# Standalone W8A8 vs BF16 MatMul benchmark (no gtest). The _biased variant
+# forces the encoding x86 uses, to show what the bias correction costs.
+add_executable(bench_matmul_i8 ops/bench_matmul_i8.cc)
+target_link_libraries(bench_matmul_i8 libgemma hwy hwy_contrib)
+
+add_executable(bench_matmul_i8_biased ops/bench_matmul_i8.cc)
+target_compile_definitions(bench_matmul_i8_biased PRIVATE
+                           GEMMA_MM_I8_FORCE_BIASED_B=1)
+target_link_libraries(bench_matmul_i8_biased libgemma hwy hwy_contrib)
+
+# W8A8 correctness, built once per A encoding so that the x86 biased-u8 path
+# is covered on non-x86 hosts as well.
+add_executable(matmul_i8_test ops/matmul_i8_test.cc)
+target_compile_definitions(matmul_i8_test PRIVATE GEMMA_MM_I8_FORCE_BIASED_B=0)
+target_link_libraries(matmul_i8_test libgemma hwy hwy_contrib)
+
+add_executable(matmul_i8_biased_test ops/matmul_i8_test.cc)
+target_compile_definitions(matmul_i8_biased_test PRIVATE
+                           GEMMA_MM_I8_FORCE_BIASED_B=1)
+target_link_libraries(matmul_i8_biased_test libgemma hwy hwy_contrib)
+
 add_executable(migrate_weights io/migrate_weights.cc)
 target_link_libraries(migrate_weights libgemma hwy hwy_contrib)
 
diff --git a/evals/compare_mmlu.py b/evals/compare_mmlu.py
new file mode 100755
index 00000000..aa7afc77
--- /dev/null
+++ b/evals/compare_mmlu.py
@@ -0,0 +1,133 @@
+#!/usr/bin/env python3
+"""Compare baseline and compressed gemma_mmlu outputs.
+
+Each input is the stdout captured from gemma_mmlu and may contain unrelated
+lines. Only lines beginning with ``MMLU_RESULT `` are parsed.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import sys
+from pathlib import Path
+from typing import Any
+
+
+RESULT_PREFIX = "MMLU_RESULT "
+
+
+def load_results(path: Path) -> dict[int, dict[str, Any]]:
+    results: dict[int, dict[str, Any]] = {}
+    with path.open(encoding="utf-8") as source:
+        for line_number, line in enumerate(source, start=1):
+            if not line.startswith(RESULT_PREFIX):
+                continue
+            try:
+                result = json.loads(line[len(RESULT_PREFIX) :])
+                question_id = int(result["id"])
+                result["correct"] = bool(result["correct"])
+                result["expected"] = str(result["expected"])
+                result["predicted"] = str(result["predicted"])
+            except (KeyError, TypeError, ValueError, json.JSONDecodeError) as error:
+                raise ValueError(
+                    f"{path}:{line_number}: invalid MMLU_RESULT: {error}"
+                ) from error
+            if question_id in results:
+                raise ValueError(f"{path}:{line_number}: duplicate id {question_id}")
+            results[question_id] = result
+
+    if not results:
+        raise ValueError(f"{path}: no {RESULT_PREFIX.strip()} lines found")
+    return results
+
+
+def compare_results(
+    baseline: dict[int, dict[str, Any]], variant: dict[int, dict[str, Any]]
+) -> dict[str, int | float]:
+    baseline_ids = set(baseline)
+    variant_ids = set(variant)
+    if baseline_ids != variant_ids:
+        missing = sorted(baseline_ids - variant_ids)
+        extra = sorted(variant_ids - baseline_ids)
+        raise ValueError(
+            "result IDs differ: "
+            f"missing from variant={missing[:10]}, extra in variant={extra[:10]}"
+        )
+
+    correct_to_incorrect = 0
+    incorrect_to_correct = 0
+    wrong_to_wrong_changes = 0
+    answer_changes = 0
+    baseline_correct = 0
+    variant_correct = 0
+
+    for question_id in sorted(baseline_ids):
+        base = baseline[question_id]
+        changed = variant[question_id]
+        if base["expected"] != changed["expected"]:
+            raise ValueError(
+                f"id {question_id}: expected answers differ: "
+                f"{base['expected']!r} != {changed['expected']!r}"
+            )
+
+        base_correct = base["correct"]
+        changed_correct = changed["correct"]
+        baseline_correct += int(base_correct)
+        variant_correct += int(changed_correct)
+        answer_changed = base["predicted"] != changed["predicted"]
+        answer_changes += int(answer_changed)
+
+        if base_correct and not changed_correct:
+            correct_to_incorrect += 1
+        elif not base_correct and changed_correct:
+            incorrect_to_correct += 1
+        elif not base_correct and not changed_correct and answer_changed:
+            wrong_to_wrong_changes += 1
+
+    samples = len(baseline_ids)
+    flips = correct_to_incorrect + incorrect_to_correct
+    return {
+        "samples": samples,
+        "baseline_correct": baseline_correct,
+        "variant_correct": variant_correct,
+        "baseline_accuracy": baseline_correct / samples,
+        "variant_accuracy": variant_correct / samples,
+        "accuracy_delta": (variant_correct - baseline_correct) / samples,
+        "correct_to_incorrect": correct_to_incorrect,
+        "incorrect_to_correct": incorrect_to_correct,
+        "flips": flips,
+        "flips_fraction": flips / samples,
+        "flips_percent": 100.0 * flips / samples,
+        "wrong_to_wrong_changes": wrong_to_wrong_changes,
+        "answer_changes": answer_changes,
+        "answer_changes_fraction": answer_changes / samples,
+        "answer_changes_percent": 100.0 * answer_changes / samples,
+    }
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(
+        description="Compare baseline and variant gemma_mmlu output streams."
+    )
+    parser.add_argument("baseline", type=Path, help="baseline gemma_mmlu stdout")
+    parser.add_argument("variant", type=Path, help="variant gemma_mmlu stdout")
+    return parser.parse_args()
+
+
+def main() -> int:
+    args = parse_args()
+    try:
+        metrics = compare_results(
+            load_results(args.baseline), load_results(args.variant)
+        )
+    except (OSError, ValueError) as error:
+        print(f"error: {error}", file=sys.stderr)
+        return 1
+
+    print(f"MMLU_FLIPS {json.dumps(metrics, sort_keys=True)}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/evals/compare_mmlu_test.py b/evals/compare_mmlu_test.py
new file mode 100755
index 00000000..6d64c1e3
--- /dev/null
+++ b/evals/compare_mmlu_test.py
@@ -0,0 +1,77 @@
+#!/usr/bin/env python3
+
+import json
+import tempfile
+import unittest
+from pathlib import Path
+
+from compare_mmlu import compare_results, load_results
+
+
+def result(question_id: int, expected: str, predicted: str) -> dict[str, object]:
+    return {
+        "id": question_id,
+        "expected": expected,
+        "predicted": predicted,
+        "correct": expected == predicted,
+    }
+
+
+class CompareMmluTest(unittest.TestCase):
+    def test_flip_counts(self) -> None:
+        baseline = {
+            1: result(1, "A", "A"),
+            2: result(2, "A", "B"),
+            3: result(3, "A", "C"),
+            4: result(4, "D", "D"),
+        }
+        variant = {
+            1: result(1, "A", "B"),
+            2: result(2, "A", "A"),
+            3: result(3, "A", "D"),
+            4: result(4, "D", "D"),
+        }
+
+        metrics = compare_results(baseline, variant)
+
+        self.assertEqual(metrics["correct_to_incorrect"], 1)
+        self.assertEqual(metrics["incorrect_to_correct"], 1)
+        self.assertEqual(metrics["flips"], 2)
+        self.assertEqual(metrics["flips_percent"], 50.0)
+        self.assertEqual(metrics["wrong_to_wrong_changes"], 1)
+        self.assertEqual(metrics["answer_changes"], 3)
+        self.assertEqual(metrics["answer_changes_percent"], 75.0)
+        self.assertEqual(metrics["accuracy_delta"], 0.0)
+
+    def test_loads_prefixed_results_and_ignores_other_lines(self) -> None:
+        with tempfile.TemporaryDirectory() as directory:
+            path = Path(directory) / "run.log"
+            rows = [result(7, "B", "B"), result(8, "C", "A")]
+            path.write_text(
+                "startup noise\n"
+                + "\n".join(f"MMLU_RESULT {json.dumps(row)}" for row in rows)
+                + "\nMMLU_SUMMARY {}\n",
+                encoding="utf-8",
+            )
+
+            loaded = load_results(path)
+
+        self.assertEqual(set(loaded), {7, 8})
+        self.assertTrue(loaded[7]["correct"])
+        self.assertFalse(loaded[8]["correct"])
+
+    def test_requires_matching_question_ids(self) -> None:
+        with self.assertRaisesRegex(ValueError, "result IDs differ"):
+            compare_results(
+                {1: result(1, "A", "A")}, {2: result(2, "A", "A")}
+            )
+
+    def test_requires_matching_expected_answers(self) -> None:
+        with self.assertRaisesRegex(ValueError, "expected answers differ"):
+            compare_results(
+                {1: result(1, "A", "A")}, {1: result(1, "B", "B")}
+            )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/evals/compare_models.py b/evals/compare_models.py
new file mode 100644
index 00000000..0efa7b36
--- /dev/null
+++ b/evals/compare_models.py
@@ -0,0 +1,361 @@
+#!/usr/bin/env python3
+"""Run generic root-vs-target model comparisons and render a report.
+
+The configuration contains arbitrary model weights, Gemma CLI arguments, and
+environment variables. No optimization (W8A8 or otherwise) is special-cased.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import os
+import re
+import subprocess
+import sys
+import time
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Any
+
+from compare_mmlu import compare_results, load_results
+
+
+@dataclass(frozen=True)
+class ModelSpec:
+    name: str
+    weights: Path
+    args: tuple[str, ...]
+    env: dict[str, str]
+
+
+@dataclass(frozen=True)
+class RunMetrics:
+    wall_seconds: float
+    peak_rss_kib: int | None
+
+
+def _resolve_path(value: str, base: Path) -> Path:
+    path = Path(value)
+    return path if path.is_absolute() else (base / path).resolve()
+
+
+def parse_model_spec(data: dict[str, Any], base: Path) -> ModelSpec:
+    if not isinstance(data, dict):
+        raise ValueError("model must be an object")
+    try:
+        name = str(data["name"])
+        weights = _resolve_path(str(data["weights"]), base)
+    except KeyError as error:
+        raise ValueError(f"model is missing {error.args[0]!r}") from error
+    if not name or re.search(r"[^A-Za-z0-9_.-]", name):
+        raise ValueError(f"invalid model name {name!r}")
+    raw_args = data.get("args", [])
+    if not isinstance(raw_args, list):
+        raise ValueError(f"{name}: args must be an array")
+    args = tuple(str(arg) for arg in raw_args)
+    raw_env = data.get("env", {})
+    if not isinstance(raw_env, dict):
+        raise ValueError(f"{name}: env must be an object")
+    env = {str(key): str(value) for key, value in raw_env.items()}
+    return ModelSpec(name=name, weights=weights, args=args, env=env)
+
+
+def _read_rss_kib(pid: int) -> int | None:
+    try:
+        status = Path(f"/proc/{pid}/status").read_text(encoding="utf-8")
+    except OSError:
+        return None
+    values: dict[str, int] = {}
+    for line in status.splitlines():
+        if line.startswith(("VmHWM:", "VmRSS:")):
+            key, value, *_ = line.split()
+            values[key.rstrip(":")] = int(value)
+    return values.get("VmHWM", values.get("VmRSS"))
+
+
+def run_command(
+    command: list[str], env_updates: dict[str, str], stdout_path: Path,
+    stderr_path: Path
+) -> RunMetrics:
+    env = os.environ.copy()
+    env.update(env_updates)
+    start = time.perf_counter()
+    peak_rss_kib: int | None = None
+    with stdout_path.open("w", encoding="utf-8") as stdout, stderr_path.open(
+        "w", encoding="utf-8"
+    ) as stderr:
+        process = subprocess.Popen(command, stdout=stdout, stderr=stderr, env=env)
+        while process.poll() is None:
+            rss = _read_rss_kib(process.pid)
+            if rss is not None:
+                peak_rss_kib = max(peak_rss_kib or 0, rss)
+            time.sleep(0.02)
+        rss = _read_rss_kib(process.pid)
+        if rss is not None:
+            peak_rss_kib = max(peak_rss_kib or 0, rss)
+        return_code = process.returncode
+    wall_seconds = time.perf_counter() - start
+    if return_code != 0:
+        tail = "\n".join(
+            stderr_path.read_text(encoding="utf-8", errors="replace").splitlines()[
+                -20:
+            ]
+        )
+        raise RuntimeError(
+            f"command failed ({return_code}): {' '.join(command)}\n{tail}"
+        )
+    return RunMetrics(wall_seconds=wall_seconds, peak_rss_kib=peak_rss_kib)
+
+
+def parse_prefixed_json(path: Path, prefix: str) -> dict[str, Any]:
+    found: dict[str, Any] | None = None
+    for line in path.read_text(encoding="utf-8").splitlines():
+        if line.startswith(prefix):
+            found = json.loads(line[len(prefix) :])
+    if found is None:
+        raise ValueError(f"{path}: no {prefix.strip()} line")
+    return found
+
+
+def parse_entropy(path: Path) -> dict[str, float | int]:
+    text = path.read_text(encoding="utf-8")
+    token_matches = re.findall(r"Number of input tokens: (\d+)", text)
+    speed_matches = re.findall(
+        r"\[([0-9.eE+-]+) tokens / sec\]", text
+    )
+    entropy_matches = re.findall(
+        r"Total cross entropy: [0-9.eE+-]+ \[cumulative: ([0-9.eE+-]+)\]",
+        text,
+    )
+    if not token_matches or not speed_matches or not entropy_matches:
+        raise ValueError(f"{path}: incomplete cross-entropy output")
+    tokens = int(token_matches[-1])
+    if tokens == 0:
+        raise ValueError(f"{path}: cross-entropy input has no tokens")
+    total_bits = float(entropy_matches[-1])
+    return {
+        "tokens": tokens,
+        "total_bits": total_bits,
+        "bits_per_token": total_bits / tokens,
+        "tokens_per_second": float(speed_matches[-1]),
+    }
+
+
+
+def render_table(rows: list[dict[str, Any]]) -> str:
+    root = rows[0]
+    root_entropy = root.get("entropy")
+    lines = [
+        "| Model | Entropy bits/token | Δ entropy | tok/s | Speedup | "
+        "MMLU accuracy | Flips | Mean KL | p95 KL | Peak RSS |",
+        "|---|---:|---:|---:|---:|---:|---:|---:|---:|---:|",
+    ]
+    for row in rows:
+        entropy = row.get("entropy")
+        if entropy and root_entropy:
+            entropy_delta = 100.0 * (
+                entropy["total_bits"] / root_entropy["total_bits"] - 1.0
+            )
+            speedup = 100.0 * (
+                entropy["tokens_per_second"]
+                / root_entropy["tokens_per_second"]
+                - 1.0
+            )
+            entropy_text = f"{entropy['bits_per_token']:.4f}"
+            delta_text = f"{entropy_delta:+.3f}%"
+            speed_text = f"{entropy['tokens_per_second']:.2f}"
+            speedup_text = f"{speedup:+.1f}%"
+        else:
+            entropy_text = delta_text = speed_text = speedup_text = "—"
+        flips = row.get("flips")
+        kl = row.get("kl")
+        rss = row.get("peak_rss_kib")
+        lines.append(
+            "| {name} | {entropy} | {delta} | {speed} | {speedup} | "
+            "{accuracy:.1f}% | {flips} | {mean_kl} | {p95_kl} | {rss} |".format(
+                name=row["name"],
+                entropy=entropy_text,
+                delta=delta_text,
+                speed=speed_text,
+                speedup=speedup_text,
+                accuracy=100.0 * row["mmlu"]["accuracy"],
+                flips="—" if flips is None else f"{flips['flips_percent']:.2f}%",
+                mean_kl="—" if kl is None else f"{kl['mean']:.6g}",
+                p95_kl="—" if kl is None else f"{kl['p95']:.6g}",
+                rss="—" if rss is None else f"{rss / 1024.0:.1f} MiB",
+            )
+        )
+    return "\n".join(lines) + "\n"
+
+
+def run_evaluation(
+    spec: ModelSpec, build_dir: Path, output_dir: Path, mmlu: Path,
+    max_questions: int, reference: Path, is_root: bool,
+    entropy_path: Path | None
+) -> tuple[dict[str, Any], Path]:
+    stem = output_dir / spec.name
+    mmlu_out = stem.with_suffix(".mmlu.out")
+    mmlu_err = stem.with_suffix(".mmlu.err")
+    command = [
+        str(build_dir / "gemma_mmlu"),
+        "--weights",
+        str(spec.weights),
+        "--input",
+        str(mmlu),
+        "--verbosity",
+        "0",
+    ]
+    if max_questions:
+        command.extend(["--max_questions", str(max_questions)])
+    command.extend(
+        ["--reference_out" if is_root else "--reference_in", str(reference)]
+    )
+    command.extend(spec.args)
+    mmlu_run = run_command(command, spec.env, mmlu_out, mmlu_err)
+    mmlu_summary = parse_prefixed_json(mmlu_out, "MMLU_SUMMARY ")
+    kl_summary = (
+        None
+        if is_root
+        else parse_prefixed_json(mmlu_out, "MMLU_KL_SUMMARY ")
+    )
+
+    entropy: dict[str, float | int] | None = None
+    entropy_run: RunMetrics | None = None
+    if entropy_path is not None:
+        entropy_out = stem.with_suffix(".entropy.out")
+        entropy_err = stem.with_suffix(".entropy.err")
+        entropy_command = [
+            str(build_dir / "single_benchmark"),
+            "--weights",
+            str(spec.weights),
+            "--cross_entropy",
+            str(entropy_path),
+            "--verbosity",
+            "0",
+            *spec.args,
+        ]
+        entropy_run = run_command(
+            entropy_command, spec.env, entropy_out, entropy_err
+        )
+        entropy = parse_entropy(entropy_out)
+
+    peak_values = [mmlu_run.peak_rss_kib]
+    if entropy_run is not None:
+        peak_values.append(entropy_run.peak_rss_kib)
+    peak_rss = max((value for value in peak_values if value is not None), default=None)
+    return (
+        {
+            "name": spec.name,
+            "weights": str(spec.weights),
+            "args": list(spec.args),
+            "env": spec.env,
+            "mmlu": mmlu_summary,
+            "kl": kl_summary,
+            "entropy": entropy,
+            "mmlu_wall_seconds": mmlu_run.wall_seconds,
+            "entropy_wall_seconds": None
+            if entropy_run is None
+            else entropy_run.wall_seconds,
+            "peak_rss_kib": peak_rss,
+        },
+        mmlu_out,
+    )
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(
+        description="Compare arbitrary target models against a root model."
+    )
+    parser.add_argument("--config", required=True, type=Path)
+    parser.add_argument("--build_dir", type=Path, default=Path("build"))
+    parser.add_argument("--output_dir", type=Path)
+    return parser.parse_args()
+
+
+def main() -> int:
+    args = parse_args()
+    try:
+        config_path = args.config.resolve()
+        config = json.loads(config_path.read_text(encoding="utf-8"))
+        base = config_path.parent
+        root = parse_model_spec(config["root"], base)
+        targets = [parse_model_spec(item, base) for item in config["targets"]]
+        if not targets:
+            raise ValueError("targets must contain at least one model")
+        names = [root.name, *(target.name for target in targets)]
+        if len(names) != len(set(names)):
+            raise ValueError("model names must be unique")
+        mmlu = _resolve_path(str(config["mmlu"]), base)
+        entropy_path = (
+            None
+            if not config.get("cross_entropy")
+            else _resolve_path(str(config["cross_entropy"]), base)
+        )
+        max_questions = int(config.get("max_questions", 0))
+        if max_questions < 0:
+            raise ValueError("max_questions must be non-negative")
+        build_dir = args.build_dir.resolve()
+        required_paths = [mmlu, root.weights]
+        required_paths.extend(target.weights for target in targets)
+        if entropy_path is not None:
+            required_paths.append(entropy_path)
+        missing = [str(path) for path in required_paths if not path.is_file()]
+        if missing:
+            raise ValueError("file does not exist: " + ", ".join(missing))
+        required_programs = [build_dir / "gemma_mmlu"]
+        if entropy_path is not None:
+            required_programs.append(build_dir / "single_benchmark")
+        missing_programs = [
+            str(path) for path in required_programs if not path.is_file()
+        ]
+        if missing_programs:
+            raise ValueError(
+                "build executable does not exist: " + ", ".join(missing_programs)
+            )
+        output_dir = (
+            args.output_dir.resolve()
+            if args.output_dir
+            else (base / f"{config_path.stem}-results").resolve()
+        )
+        output_dir.mkdir(parents=True, exist_ok=True)
+        reference = output_dir / f"{root.name}.root-kl.bin"
+
+        root_row, root_output = run_evaluation(
+            root, build_dir, output_dir, mmlu, max_questions,
+            reference, True, entropy_path
+        )
+        root_row["flips"] = None
+        root_results = load_results(root_output)
+        rows = [root_row]
+        for target in targets:
+            row, target_output = run_evaluation(
+                target, build_dir, output_dir, mmlu,
+                max_questions, reference, False, entropy_path
+            )
+            row["flips"] = compare_results(
+                root_results, load_results(target_output)
+            )
+            rows.append(row)
+
+        report = {
+            "schema_version": 1,
+            "mmlu": str(mmlu),
+            "cross_entropy": None if entropy_path is None else str(entropy_path),
+            "reference": str(reference),
+            "models": rows,
+        }
+        (output_dir / "comparison.json").write_text(
+            json.dumps(report, indent=2, sort_keys=True) + "\n", encoding="utf-8"
+        )
+        table = render_table(rows)
+        (output_dir / "comparison.md").write_text(table, encoding="utf-8")
+        print(table, end="")
+        return 0
+    except (KeyError, OSError, ValueError, RuntimeError, json.JSONDecodeError) as error:
+        print(f"error: {error}", file=sys.stderr)
+        return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/evals/compare_models_test.py b/evals/compare_models_test.py
new file mode 100644
index 00000000..112e03e1
--- /dev/null
+++ b/evals/compare_models_test.py
@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+
+import tempfile
+import unittest
+from pathlib import Path
+
+from compare_models import (
+    parse_entropy,
+    parse_model_spec,
+    parse_prefixed_json,
+    render_table,
+)
+
+
+class CompareModelsTest(unittest.TestCase):
+    def test_parse_model_spec_resolves_paths_and_environment(self) -> None:
+        spec = parse_model_spec(
+            {
+                "name": "w8a8",
+                "weights": "models/model.sbs",
+                "args": ["--num_threads", 4],
+                "env": {"GEMMA_MM_I8": 1},
+            },
+            Path("/work"),
+        )
+
+        self.assertEqual(spec.name, "w8a8")
+        self.assertEqual(spec.weights, Path("/work/models/model.sbs"))
+        self.assertEqual(spec.args, ("--num_threads", "4"))
+        self.assertEqual(spec.env, {"GEMMA_MM_I8": "1"})
+
+    def test_parse_prefixed_json_uses_last_summary(self) -> None:
+        with tempfile.TemporaryDirectory() as directory:
+            path = Path(directory) / "output.log"
+            path.write_text(
+                'MMLU_SUMMARY {"answers":1}\n'
+                'noise\nMMLU_SUMMARY {"answers":2}\n',
+                encoding="utf-8",
+            )
+            parsed = parse_prefixed_json(path, "MMLU_SUMMARY ")
+
+        self.assertEqual(parsed, {"answers": 2})
+
+    def test_parse_entropy(self) -> None:
+        with tempfile.TemporaryDirectory() as directory:
+            path = Path(directory) / "entropy.log"
+            path.write_text(
+                "Number of input tokens: 8\n"
+                "Took 1.0 s [8.0 tokens / sec]\n"
+                "Total cross entropy: 12.0 [cumulative: 12.0]\n",
+                encoding="utf-8",
+            )
+            parsed = parse_entropy(path)
+
+        self.assertEqual(parsed["tokens"], 8)
+        self.assertEqual(parsed["total_bits"], 12.0)
+        self.assertEqual(parsed["bits_per_token"], 1.5)
+        self.assertEqual(parsed["tokens_per_second"], 8.0)
+
+    def test_render_table(self) -> None:
+        root = {
+            "name": "root",
+            "mmlu": {"accuracy": 0.5},
+            "entropy": {
+                "total_bits": 20.0,
+                "bits_per_token": 2.0,
+                "tokens_per_second": 10.0,
+            },
+            "flips": None,
+            "kl": None,
+            "peak_rss_kib": 1024,
+        }
+        target = {
+            "name": "target",
+            "mmlu": {"accuracy": 0.75},
+            "entropy": {
+                "total_bits": 22.0,
+                "bits_per_token": 2.2,
+                "tokens_per_second": 12.0,
+            },
+            "flips": {"flips_percent": 25.0},
+            "kl": {"mean": 0.01, "p95": 0.03},
+            "peak_rss_kib": 2048,
+        }
+
+        table = render_table([root, target])
+
+        self.assertIn("| root | 2.0000 | +0.000% | 10.00 | +0.0%", table)
+        self.assertIn(
+            "| target | 2.2000 | +10.000% | 12.00 | +20.0% | "
+            "75.0% | 25.00% | 0.01 | 0.03 | 2.0 MiB |",
+            table,
+        )
+
+    def test_rejects_unsafe_report_name(self) -> None:
+        with self.assertRaisesRegex(ValueError, "invalid model name"):
+            parse_model_spec(
+                {"name": "../target", "weights": "model.sbs"}, Path("/work")
+            )
+
+    def test_rejects_string_model_args(self) -> None:
+        with self.assertRaisesRegex(ValueError, "args must be an array"):
+            parse_model_spec(
+                {"name": "target", "weights": "model.sbs", "args": "--foo"},
+                Path("/work"),
+            )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/evals/model_comparison.cc b/evals/model_comparison.cc
new file mode 100644
index 00000000..cb61c30d
--- /dev/null
+++ b/evals/model_comparison.cc
@@ -0,0 +1,189 @@
+// Copyright 2026 Google LLC
+// SPDX-License-Identifier: Apache-2.0
+
+#include "evals/model_comparison.h"
+
+#include <algorithm>
+#include <cmath>
+#include <cstring>
+#include <limits>
+#include <stdexcept>
+#include <type_traits>
+
+namespace gcpp {
+namespace {
+
+constexpr char kMagic[8] = {'G', 'C', 'P', 'P', 'K', 'L', '0', '1'};
+constexpr uint32_t kVersion = 1;
+constexpr uint32_t kEndianMarker = 0x01020304u;
+
+template <typename T>
+void WriteValue(std::ofstream& stream, const T& value) {
+  static_assert(std::is_trivially_copyable<T>::value, "binary scalar");
+  stream.write(reinterpret_cast<const char*>(&value), sizeof(value));
+  if (!stream) throw std::runtime_error("failed to write KL reference file");
+}
+
+template <typename T>
+T ReadValue(std::ifstream& stream) {
+  static_assert(std::is_trivially_copyable<T>::value, "binary scalar");
+  T value;
+  stream.read(reinterpret_cast<char*>(&value), sizeof(value));
+  if (!stream) throw std::runtime_error("truncated KL reference file");
+  return value;
+}
+
+void RequireEqual(const char* name, uint64_t actual, uint64_t expected) {
+  if (actual != expected) {
+    throw std::runtime_error(std::string("KL reference ") + name +
+                             " mismatch: " + std::to_string(actual) +
+                             " != " + std::to_string(expected));
+  }
+}
+
+}  // namespace
+
+uint64_t ModelComparisonFingerprint(const std::string& bytes) {
+  uint64_t hash = 14695981039346656037ull;
+  for (const unsigned char byte : bytes) {
+    hash ^= byte;
+    hash *= 1099511628211ull;
+  }
+  return hash;
+}
+
+double FullVocabLogSumExp(const float* logits, size_t size) {
+  if (size == 0) throw std::invalid_argument("empty logits");
+  float max_logit = -std::numeric_limits<float>::infinity();
+  for (size_t i = 0; i < size; ++i) {
+    max_logit = std::max(max_logit, logits[i]);
+  }
+  if (!std::isfinite(max_logit)) {
+    throw std::invalid_argument("non-finite maximum logit");
+  }
+
+  double sum = 0.0;
+  for (size_t i = 0; i < size; ++i) {
+    sum += std::exp(static_cast<double>(logits[i] - max_logit));
+  }
+  return static_cast<double>(max_logit) + std::log(sum);
+}
+
+double FullVocabKLDivergence(const std::vector<float>& root_logits,
+                             double root_log_sum_exp,
+                             const float* target_logits, size_t size) {
+  if (root_logits.size() != size) {
+    throw std::invalid_argument("root/target vocabulary size mismatch");
+  }
+  const double target_log_sum_exp = FullVocabLogSumExp(target_logits, size);
+  double kl = 0.0;
+  for (size_t i = 0; i < size; ++i) {
+    const double root_log_prob =
+        static_cast<double>(root_logits[i]) - root_log_sum_exp;
+    const double target_log_prob =
+        static_cast<double>(target_logits[i]) - target_log_sum_exp;
+    kl += std::exp(root_log_prob) * (root_log_prob - target_log_prob);
+  }
+  return kl < 0.0 && kl > -1E-12 ? 0.0 : kl;
+}
+
+ModelComparisonWriter::ModelComparisonWriter(
+    const std::string& path, const ModelComparisonMetadata& metadata)
+    : stream_(path, std::ios::binary), metadata_(metadata) {
+  static_assert(sizeof(float) == 4, "reference format requires float32");
+  if (!stream_) throw std::runtime_error("cannot create KL reference: " + path);
+  stream_.write(kMagic, sizeof(kMagic));
+  WriteValue(stream_, kVersion);
+  WriteValue(stream_, kEndianMarker);
+  WriteValue(stream_, metadata_.vocab_size);
+  WriteValue(stream_, metadata_.sample_count);
+  WriteValue(stream_, metadata_.dataset_fingerprint);
+  WriteValue(stream_, metadata_.tokenizer_fingerprint);
+}
+
+ModelComparisonWriter::~ModelComparisonWriter() {
+  if (!finished_) stream_.close();
+}
+
+void ModelComparisonWriter::Write(int64_t sample_id, int32_t expected_label,
+                                  const float* logits, size_t size) {
+  if (finished_) throw std::runtime_error("KL reference already finished");
+  if (size != metadata_.vocab_size) {
+    throw std::invalid_argument("logits do not match reference vocabulary");
+  }
+  if (records_written_ >= metadata_.sample_count) {
+    throw std::runtime_error("too many KL reference records");
+  }
+  WriteValue(stream_, sample_id);
+  WriteValue(stream_, expected_label);
+  const double log_sum_exp = FullVocabLogSumExp(logits, size);
+  WriteValue(stream_, log_sum_exp);
+  stream_.write(reinterpret_cast<const char*>(logits),
+                static_cast<std::streamsize>(size * sizeof(float)));
+  if (!stream_) throw std::runtime_error("failed to write KL reference logits");
+  ++records_written_;
+}
+
+void ModelComparisonWriter::Finish() {
+  if (finished_) return;
+  if (records_written_ != metadata_.sample_count) {
+    throw std::runtime_error("KL reference record count mismatch");
+  }
+  stream_.flush();
+  if (!stream_) throw std::runtime_error("failed to finish KL reference file");
+  finished_ = true;
+}
+
+ModelComparisonReader::ModelComparisonReader(const std::string& path)
+    : stream_(path, std::ios::binary) {
+  if (!stream_) throw std::runtime_error("cannot open KL reference: " + path);
+  char magic[sizeof(kMagic)];
+  stream_.read(magic, sizeof(magic));
+  if (!stream_ || std::memcmp(magic, kMagic, sizeof(kMagic)) != 0) {
+    throw std::runtime_error("invalid KL reference magic");
+  }
+  RequireEqual("version", ReadValue<uint32_t>(stream_), kVersion);
+  RequireEqual("endianness", ReadValue<uint32_t>(stream_), kEndianMarker);
+  metadata_.vocab_size = ReadValue<uint32_t>(stream_);
+  metadata_.sample_count = ReadValue<uint64_t>(stream_);
+  metadata_.dataset_fingerprint = ReadValue<uint64_t>(stream_);
+  metadata_.tokenizer_fingerprint = ReadValue<uint64_t>(stream_);
+}
+
+void ModelComparisonReader::Validate(
+    const ModelComparisonMetadata& expected) const {
+  RequireEqual("vocabulary", metadata_.vocab_size, expected.vocab_size);
+  RequireEqual("sample count", metadata_.sample_count, expected.sample_count);
+  RequireEqual("dataset fingerprint", metadata_.dataset_fingerprint,
+               expected.dataset_fingerprint);
+  RequireEqual("tokenizer fingerprint", metadata_.tokenizer_fingerprint,
+               expected.tokenizer_fingerprint);
+}
+
+ModelComparisonRecord ModelComparisonReader::Read() {
+  if (records_read_ >= metadata_.sample_count) {
+    throw std::runtime_error("too many KL reference reads");
+  }
+  ModelComparisonRecord record;
+  record.sample_id = ReadValue<int64_t>(stream_);
+  record.expected_label = ReadValue<int32_t>(stream_);
+  record.log_sum_exp = ReadValue<double>(stream_);
+  record.logits.resize(metadata_.vocab_size);
+  stream_.read(
+      reinterpret_cast<char*>(record.logits.data()),
+      static_cast<std::streamsize>(record.logits.size() * sizeof(float)));
+  if (!stream_) throw std::runtime_error("truncated KL reference logits");
+  ++records_read_;
+  return record;
+}
+
+void ModelComparisonReader::Finish() {
+  if (records_read_ != metadata_.sample_count) {
+    throw std::runtime_error("unread KL reference records");
+  }
+  if (stream_.peek() != std::ifstream::traits_type::eof()) {
+    throw std::runtime_error("trailing bytes in KL reference file");
+  }
+}
+
+}  // namespace gcpp
diff --git a/evals/model_comparison.h b/evals/model_comparison.h
new file mode 100644
index 00000000..cb8b276f
--- /dev/null
+++ b/evals/model_comparison.h
@@ -0,0 +1,75 @@
+// Copyright 2026 Google LLC
+// SPDX-License-Identifier: Apache-2.0
+
+#ifndef THIRD_PARTY_GEMMA_CPP_EVALS_MODEL_COMPARISON_H_
+#define THIRD_PARTY_GEMMA_CPP_EVALS_MODEL_COMPARISON_H_
+
+#include <stdint.h>
+
+#include <cstddef>
+#include <fstream>
+#include <string>
+#include <vector>
+
+namespace gcpp {
+
+uint64_t ModelComparisonFingerprint(const std::string& bytes);
+
+// Numerically stable full-vocabulary operations. KL is directional:
+// D_KL(root || target).
+double FullVocabLogSumExp(const float* logits, size_t size);
+double FullVocabKLDivergence(const std::vector<float>& root_logits,
+                             double root_log_sum_exp,
+                             const float* target_logits, size_t size);
+
+struct ModelComparisonMetadata {
+  uint32_t vocab_size = 0;
+  uint64_t sample_count = 0;
+  uint64_t dataset_fingerprint = 0;
+  uint64_t tokenizer_fingerprint = 0;
+};
+
+struct ModelComparisonRecord {
+  int64_t sample_id = 0;
+  int32_t expected_label = 0;
+  double log_sum_exp = 0.0;
+  std::vector<float> logits;
+};
+
+// Versioned binary store for root-model logits. It is uncompressed so target
+// runs can stream one question at a time without loading the whole dataset.
+class ModelComparisonWriter {
+ public:
+  ModelComparisonWriter(const std::string& path,
+                        const ModelComparisonMetadata& metadata);
+  ~ModelComparisonWriter();
+
+  void Write(int64_t sample_id, int32_t expected_label, const float* logits,
+             size_t size);
+  void Finish();
+
+ private:
+  std::ofstream stream_;
+  ModelComparisonMetadata metadata_;
+  uint64_t records_written_ = 0;
+  bool finished_ = false;
+};
+
+class ModelComparisonReader {
+ public:
+  explicit ModelComparisonReader(const std::string& path);
+
+  const ModelComparisonMetadata& Metadata() const { return metadata_; }
+  void Validate(const ModelComparisonMetadata& expected) const;
+  ModelComparisonRecord Read();
+  void Finish();
+
+ private:
+  std::ifstream stream_;
+  ModelComparisonMetadata metadata_;
+  uint64_t records_read_ = 0;
+};
+
+}  // namespace gcpp
+
+#endif  // THIRD_PARTY_GEMMA_CPP_EVALS_MODEL_COMPARISON_H_
diff --git a/evals/model_comparison_test.cc b/evals/model_comparison_test.cc
new file mode 100644
index 00000000..57ed1962
--- /dev/null
+++ b/evals/model_comparison_test.cc
@@ -0,0 +1,86 @@
+// Copyright 2026 Google LLC
+// SPDX-License-Identifier: Apache-2.0
+
+#include "evals/model_comparison.h"
+
+#include <stdio.h>
+
+#include <cmath>
+#include <cstdio>
+#include <vector>
+
+namespace {
+
+int failures = 0;
+
+void CheckNear(const char* name, double actual, double expected,
+               double tolerance = 1E-12) {
+  if (std::abs(actual - expected) > tolerance) {
+    fprintf(stderr, "FAIL %s: %.17g != %.17g\n", name, actual, expected);
+    ++failures;
+  }
+}
+
+void TestKL() {
+  const std::vector<float> root = {std::log(0.25f), std::log(0.75f)};
+  const std::vector<float> same = root;
+  const std::vector<float> shifted = {root[0] + 17.0f, root[1] + 17.0f};
+  const std::vector<float> uniform = {0.0f, 0.0f};
+  const double root_lse = gcpp::FullVocabLogSumExp(root.data(), root.size());
+
+  CheckNear(
+      "identical",
+      gcpp::FullVocabKLDivergence(root, root_lse, same.data(), same.size()),
+      0.0);
+  CheckNear("shift invariant",
+            gcpp::FullVocabKLDivergence(root, root_lse, shifted.data(),
+                                        shifted.size()),
+            0.0, 1E-7);
+  const double expected = 0.25 * std::log(0.5) + 0.75 * std::log(1.5);
+  CheckNear("known KL",
+            gcpp::FullVocabKLDivergence(root, root_lse, uniform.data(),
+                                        uniform.size()),
+            expected, 1E-7);
+}
+
+void TestReferenceRoundTrip() {
+  const char* path = "/tmp/gemma_model_comparison_test.bin";
+  std::remove(path);
+  const gcpp::ModelComparisonMetadata metadata = {
+      /*vocab_size=*/3,
+      /*sample_count=*/1,
+      /*dataset_fingerprint=*/123,
+      /*tokenizer_fingerprint=*/456,
+  };
+  const std::vector<float> logits = {1.0f, -2.0f, 4.0f};
+  {
+    gcpp::ModelComparisonWriter writer(path, metadata);
+    writer.Write(7, 2, logits.data(), logits.size());
+    writer.Finish();
+  }
+  {
+    gcpp::ModelComparisonReader reader(path);
+    reader.Validate(metadata);
+    const gcpp::ModelComparisonRecord record = reader.Read();
+    if (record.sample_id != 7 || record.expected_label != 2 ||
+        record.logits != logits) {
+      fprintf(stderr, "FAIL reference round trip\n");
+      ++failures;
+    }
+    reader.Finish();
+  }
+  std::remove(path);
+}
+
+}  // namespace
+
+int main() {
+  TestKL();
+  TestReferenceRoundTrip();
+  if (failures != 0) {
+    fprintf(stderr, "FAIL (%d failures)\n", failures);
+    return 1;
+  }
+  printf("PASS\n");
+  return 0;
+}
diff --git a/evals/run_mmlu.cc b/evals/run_mmlu.cc
index 66044397..0b448153 100644
--- a/evals/run_mmlu.cc
+++ b/evals/run_mmlu.cc
@@ -16,17 +16,24 @@
 #include <stdio.h>
 
 #include <algorithm>
+#include <array>
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <numeric>
+#include <stdexcept>
 #include <string>
 #include <vector>
 
 #include "evals/benchmark_helper.h"
+#include "evals/model_comparison.h"
 #include "gemma/gemma.h"  // Gemma
-#include "io/io.h"        // Path
-#include "util/args.h"
 #include "hwy/base.h"
 #include "hwy/highway.h"
 #include "hwy/profiler.h"
+#include "io/io.h"  // Path
 #include "nlohmann/json.hpp"
+#include "util/args.h"
 
 namespace gcpp {
 
@@ -36,70 +43,125 @@ struct JsonArgs : public ArgsBase<JsonArgs> {
   }
 
   Path input;
+  Path reference_out;
+  Path reference_in;
+  size_t max_questions;
 
-  // Returns error string or nullptr if OK.
   const char* Validate() const {
     if (input.Empty()) return "Must specify --input";
     if (!input.Exists()) return "--input file does not exist";
+    if (!reference_out.Empty() && !reference_in.Empty()) {
+      return "Specify only one of --reference_out and --reference_in";
+    }
+    if (!reference_in.Empty() && !reference_in.Exists()) {
+      return "--reference_in file does not exist";
+    }
     return nullptr;
   }
 
   template <class Visitor>
   void ForEach(const Visitor& visitor) {
     visitor(input, "input", Path(), "Full pathname of mmlu.json.");
-  };
+    visitor(reference_out, "reference_out", Path(),
+            "Write root-model full-vocabulary logits to this binary file.");
+    visitor(reference_in, "reference_in", Path(),
+            "Compare this target model against a root reference file.");
+    visitor(max_questions, "max_questions", size_t{0},
+            "Maximum questions to run; zero runs the full dataset.");
+  }
 };
 
-// Linear search for a few tokens is faster than std::set.
-// TODO: instead of accepting for each vocab entry, filter the logits once.
-class TokenSet {
+// Maps both "A" and " A" tokenizer variants to answer labels 0..3.
+class AnswerTokens {
  public:
-  TokenSet(const GemmaTokenizer& tokenizer,
-           const std::vector<std::string>& strings) {
-    all_tokens_.reserve(strings.size());
-    for (const std::string& str : strings) {
-      std::vector<int> tokens;
-      fprintf(stderr, "%s -> ", str.c_str());
-      HWY_ASSERT(tokenizer.Encode(str, &tokens));
-      for (int token : tokens) {
-        fprintf(stderr, "%d, ", token);
-        all_tokens_.push_back(token);
+  explicit AnswerTokens(const GemmaTokenizer& tokenizer) {
+    for (int label = 0; label < 4; ++label) {
+      for (const std::string& prefix : {std::string(), std::string(" ")}) {
+        const std::string str = prefix + static_cast<char>('A' + label);
+        std::vector<int> tokens;
+        HWY_ASSERT(tokenizer.Encode(str, &tokens));
+        HWY_ASSERT(tokens.size() == 1);
+        fprintf(stderr, "%s -> %d\n", str.c_str(), tokens[0]);
+        tokens_.push_back({tokens[0], label});
       }
-      fprintf(stderr, "\n");
     }
   }
 
-  bool Contains(int token) const {
-    return std::find(all_tokens_.begin(), all_tokens_.end(), token) !=
-           all_tokens_.end();
+  int Label(int token) const {
+    const auto it =
+        std::find_if(tokens_.begin(), tokens_.end(),
+                     [token](const auto& item) { return item.first == token; });
+    return it == tokens_.end() ? -1 : it->second;
   }
 
+  const std::vector<std::pair<int, int>>& All() const { return tokens_; }
+
  private:
-  std::vector<int> all_tokens_;
+  std::vector<std::pair<int, int>> tokens_;
 };
 
-void Run(GemmaEnv& env, JsonArgs& json) {
+double Percentile(const std::vector<double>& sorted, double quantile) {
+  if (sorted.empty()) return 0.0;
+  const double position = quantile * static_cast<double>(sorted.size() - 1);
+  const size_t lower = static_cast<size_t>(std::floor(position));
+  const size_t upper = static_cast<size_t>(std::ceil(position));
+  const double fraction = position - static_cast<double>(lower);
+  return sorted[lower] + fraction * (sorted[upper] - sorted[lower]);
+}
+
+void Run(GemmaEnv& env, JsonArgs& args) {
   PROFILER_ZONE("Run.all");
 
-  float answers = 0.0f;
-  float correct_answers = 0.0f;
+  size_t answers = 0;
+  size_t correct_answers = 0;
+  std::vector<double> kl_values;
+
+  const std::string json_text = ReadFileToString(args.input);
+  const auto json_data = nlohmann::json::parse(json_text);
+  const auto& samples = json_data["samples"];
+  const size_t sample_count =
+      args.max_questions == 0
+          ? samples.size()
+          : std::min<size_t>(args.max_questions, samples.size());
+
+  const Gemma& gemma = *env.GetGemma();
+  const GemmaTokenizer& tokenizer = gemma.Tokenizer();
+  const ModelComparisonMetadata metadata = {
+      static_cast<uint32_t>(gemma.Config().vocab_size), sample_count,
+      ModelComparisonFingerprint(json_text),
+      ModelComparisonFingerprint(tokenizer.Serialize())};
+
+  std::unique_ptr<ModelComparisonWriter> reference_writer;
+  std::unique_ptr<ModelComparisonReader> reference_reader;
+  if (!args.reference_out.Empty()) {
+    reference_writer = std::make_unique<ModelComparisonWriter>(
+        args.reference_out.path, metadata);
+  } else if (!args.reference_in.Empty()) {
+    reference_reader =
+        std::make_unique<ModelComparisonReader>(args.reference_in.path);
+    reference_reader->Validate(metadata);
+    kl_values.reserve(sample_count);
+  }
+
+  const AnswerTokens answer_tokens(tokenizer);
 
-  auto json_data = nlohmann::json::parse(ReadFileToString(json.input));
+  for (const auto& sample : samples) {
+    if (answers >= sample_count) break;
+    const int64_t id = sample["i"];
+    fprintf(stderr, "Processing question %lld\n", static_cast<long long>(id));
+    const int correct_label = sample["input_label"];
+    const std::string correct_answer(1, static_cast<char>('A' + correct_label));
 
-  const std::vector<std::string> accept_strings = {
-      "A",  "B",   "C",   "D",   //
-      " A", " B",  " C",  " D",  //
-      "**", "**:", ":**", "The", "Answer", "is", ":", "."};
-  const TokenSet accept_set(env.GetGemma()->Tokenizer(), accept_strings);
+    ModelComparisonRecord root_record;
+    if (reference_reader) {
+      root_record = reference_reader->Read();
+      if (root_record.sample_id != id ||
+          root_record.expected_label != correct_label) {
+        throw std::runtime_error("KL reference sample identity mismatch");
+      }
+    }
 
-  for (auto sample : json_data["samples"]) {
-    const int id = sample["i"];
-    fprintf(stderr, "Processing question %d\n", id);
-    const std::string& correct_answer = accept_strings[sample["input_label"]];
     std::string prompt_string = sample["prompt"];
-    // AcceptFunc restricts the output to one of these four tokens, so make an
-    // effort to steer the model towards that. See
-    // https://huggingface.co/blog/open-llm-leaderboard-mmlu
     prompt_string +=
         "What is start of the line with the correct answer? "
         "Do not include any justifications or explanations. Reply only with a "
@@ -107,45 +169,149 @@ void Run(GemmaEnv& env, JsonArgs& json) {
     const std::vector<int> prompt = env.WrapAndTokenize(prompt_string);
     const size_t prompt_size = prompt.size();
 
-    std::vector<int> predicted_token_ids;
-    predicted_token_ids.reserve(4096);
+    int predicted_token = -1;
+    std::array<float, 4> answer_logits;
+    std::array<float, 4> answer_probs;
+    answer_logits.fill(-std::numeric_limits<float>::infinity());
+    answer_probs.fill(0.0f);
+    std::vector<float> captured_logits;
+    double full_vocab_kl = 0.0;
     size_t generated = 0;
-    const StreamFunc stream_token = [&generated, prompt_size,
-                                     &predicted_token_ids](int token,
-                                                           float proba) {
+    const StreamFunc stream_token = [&generated, prompt_size, &predicted_token](
+                                        int token, float /*proba*/) {
       PROFILER_ZONE("Stream");
       ++generated;
       if (generated > prompt_size) {
-        predicted_token_ids.push_back(token);
+        predicted_token = token;
+        return false;
       }
       return true;
     };
 
-    // Although " A" is a token, it is difficult to associate that with the
-    // correct answer. Only accepting certain tokens is risky: (A) is easily
-    // confused with the word "A".
     gcpp::TimingInfo timing_info;
     gcpp::RuntimeConfig runtime_config = {
-        .max_generated_tokens = 30,
+        .max_generated_tokens = 1,
         .temperature = 0.0f,
         .verbosity = env.Verbosity(),
         .attention_impl = env.MutableConfig().attention_impl,
         .stream_token = stream_token,
+        .sample_func = [&answer_tokens, &answer_logits, &answer_probs,
+                        &captured_logits, &full_vocab_kl, &reference_writer,
+                        &reference_reader, &root_record](
+                           size_t /*query_idx*/, size_t /*pos*/, Logits logits,
+                           size_t /*worker*/) -> TokenAndProb {
+          if (reference_writer) {
+            captured_logits.assign(logits.data(),
+                                   logits.data() + logits.size());
+          } else if (reference_reader) {
+            full_vocab_kl = FullVocabKLDivergence(root_record.logits,
+                                                  root_record.log_sum_exp,
+                                                  logits.data(), logits.size());
+          }
+
+          int best_token = -1;
+          int best_label = -1;
+          float best_logit = -std::numeric_limits<float>::infinity();
+          for (const auto& [token, label] : answer_tokens.All()) {
+            if (logits[token] > answer_logits[label]) {
+              answer_logits[label] = logits[token];
+            }
+            if (logits[token] > best_logit) {
+              best_logit = logits[token];
+              best_token = token;
+              best_label = label;
+            }
+          }
+
+          float sum = 0.0f;
+          for (int label = 0; label < 4; ++label) {
+            answer_probs[label] = std::exp(answer_logits[label] - best_logit);
+            sum += answer_probs[label];
+          }
+          for (float& prob : answer_probs) prob /= sum;
+          return TokenAndProb{.token = best_token,
+                              .prob = answer_probs[best_label]};
+        },
     };
     env.GetGemma()->Generate(runtime_config, prompt, /*pos=*/0,
                              env.MutableKVCache(), env.MutableEnv(),
                              timing_info);
 
-    std::string output_string = env.StringFromTokens(predicted_token_ids);
+    if (reference_writer) {
+      if (captured_logits.size() != metadata.vocab_size) {
+        throw std::runtime_error(
+            "failed to capture root full-vocabulary logits");
+      }
+      reference_writer->Write(id, correct_label, captured_logits.data(),
+                              captured_logits.size());
+    } else if (reference_reader) {
+      if (!std::isfinite(full_vocab_kl) || full_vocab_kl < 0.0) {
+        throw std::runtime_error("invalid full-vocabulary KL divergence");
+      }
+      kl_values.push_back(full_vocab_kl);
+    }
+
+    const int predicted_label = answer_tokens.Label(predicted_token);
+    const std::string output_string =
+        predicted_label == -1
+            ? std::string("?")
+            : std::string(1, static_cast<char>('A' + predicted_label));
     fprintf(stderr, "Correct %s, model '%s'\n", correct_answer.c_str(),
             output_string.c_str());
 
-    answers += 1.0f;
-    if (output_string == correct_answer) {
-      correct_answers += 1.0f;
+    const bool is_correct = predicted_label == correct_label;
+    float second_logit = -std::numeric_limits<float>::infinity();
+    for (int label = 0; label < 4; ++label) {
+      if (label != predicted_label) {
+        second_logit = std::max(second_logit, answer_logits[label]);
+      }
     }
-    fprintf(stderr, "%.0f/%.0f = %.2f%%\n", correct_answers, answers,
-            correct_answers / answers);
+    ++answers;
+    correct_answers += static_cast<size_t>(is_correct);
+    nlohmann::json result = {
+        {"id", id},
+        {"expected", correct_answer},
+        {"predicted", output_string},
+        {"correct", is_correct},
+        {"logits", answer_logits},
+        {"probabilities", answer_probs},
+        {"margin", predicted_label == -1
+                       ? 0.0f
+                       : answer_logits[predicted_label] - second_logit},
+    };
+    if (reference_reader) result["full_vocab_kl"] = full_vocab_kl;
+    printf("MMLU_RESULT %s\n", result.dump().c_str());
+    fflush(stdout);
+    fprintf(stderr, "%zu/%zu = %.2f%%\n", correct_answers, answers,
+            100.0 * static_cast<double>(correct_answers) / answers);
+  }
+
+  if (reference_writer) reference_writer->Finish();
+  if (reference_reader) reference_reader->Finish();
+
+  const nlohmann::json summary = {
+      {"answers", answers},
+      {"correct", correct_answers},
+      {"accuracy",
+       answers == 0 ? 0.0 : static_cast<double>(correct_answers) / answers},
+  };
+  printf("MMLU_SUMMARY %s\n", summary.dump().c_str());
+
+  if (!kl_values.empty()) {
+    std::sort(kl_values.begin(), kl_values.end());
+    const double mean =
+        std::accumulate(kl_values.begin(), kl_values.end(), 0.0) /
+        kl_values.size();
+    const nlohmann::json kl_summary = {
+        {"samples", kl_values.size()},
+        {"mean", mean},
+        {"median", Percentile(kl_values, 0.5)},
+        {"p95", Percentile(kl_values, 0.95)},
+        {"max", kl_values.back()},
+        {"unit", "nats"},
+        {"direction", "root||target"},
+    };
+    printf("MMLU_KL_SUMMARY %s\n", kl_summary.dump().c_str());
   }
 }
 
@@ -154,17 +320,22 @@ void Run(GemmaEnv& env, JsonArgs& json) {
 int main(int argc, char** argv) {
   gcpp::InternalInit();
 
-  {
-    PROFILER_ZONE("Startup.all");
-    gcpp::ConsumedArgs consumed(argc, argv);
-    gcpp::GemmaArgs args(argc, argv, consumed);
-    gcpp::JsonArgs json_args(argc, argv, consumed);
-    gcpp::AbortIfInvalidArgs(json_args);
-    consumed.AbortIfUnconsumed();
+  try {
+    {
+      PROFILER_ZONE("Startup.all");
+      gcpp::ConsumedArgs consumed(argc, argv);
+      gcpp::GemmaArgs args(argc, argv, consumed);
+      gcpp::JsonArgs json_args(argc, argv, consumed);
+      gcpp::AbortIfInvalidArgs(json_args);
+      consumed.AbortIfUnconsumed();
 
-    gcpp::GemmaEnv env(args);
-    gcpp::Run(env, json_args);
+      gcpp::GemmaEnv env(args);
+      gcpp::Run(env, json_args);
+    }
+    PROFILER_PRINT_RESULTS();
+    return 0;
+  } catch (const std::exception& error) {
+    fprintf(stderr, "model comparison failed: %s\n", error.what());
+    return 1;
   }
-  PROFILER_PRINT_RESULTS();  // Must call outside the zone above.
-  return 0;
 }
diff --git a/ops/bench_matmul_i8.cc b/ops/bench_matmul_i8.cc
new file mode 100644
index 00000000..825db78e
--- /dev/null
+++ b/ops/bench_matmul_i8.cc
@@ -0,0 +1,399 @@
+// Copyright 2025 Google LLC
+// SPDX-License-Identifier: Apache-2.0
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//      http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Compares the BF16 MatMul (`ops/matmul-inl.h`) against the W8A8 int8 kernel
+// (`ops/matmul_i8-inl.h`) on Gemma-shaped problems, and reports the accuracy
+// of both relative to an F64 reference. Standalone binary (no gtest) so that
+// it can be run directly.
+
+#include <math.h>
+
+#include <cmath>
+#include <stddef.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <string.h>
+
+#include <algorithm>
+#include <vector>
+
+#include "compression/types.h"  // GEMMA_DISABLED_TARGETS
+#ifndef HWY_DISABLED_TARGETS
+#define HWY_DISABLED_TARGETS GEMMA_DISABLED_TARGETS
+#endif  // HWY_DISABLED_TARGETS
+
+#include "ops/matmul.h"
+#include "util/basics.h"
+#include "util/mat.h"
+#include "util/threading_context.h"
+#include "hwy/aligned_allocator.h"
+#include "hwy/timer.h"
+
+// clang-format off
+#undef HWY_TARGET_INCLUDE
+#define HWY_TARGET_INCLUDE "ops/bench_matmul_i8.cc"  // NOLINT
+// clang-format on
+#include "hwy/foreach_target.h"  // IWYU pragma: keep
+#include "hwy/highway.h"
+// After highway.h
+#include "compression/compress-inl.h"
+#include "ops/matmul-inl.h"
+#include "ops/matmul_i8-inl.h"
+
+HWY_BEFORE_NAMESPACE();
+namespace gcpp {
+namespace HWY_NAMESPACE {
+namespace hn = hwy::HWY_NAMESPACE;
+
+// Deterministic, reproducible pseudo-Gaussian. Real activations and weights
+// are roughly bell-shaped; the ramp in `compression/test_util-inl.h` would
+// flatter or penalize int8 quantization for the wrong reasons.
+class Rng {
+ public:
+  explicit Rng(uint64_t seed) : state_(seed * 6364136223846793005ull + 1) {}
+
+  float Normal() {
+    // Sum of 4 uniforms: close enough to Gaussian, and cheap.
+    float sum = 0.0f;
+    for (int i = 0; i < 4; ++i) sum += Uniform();
+    return (sum - 2.0f) * 1.732f;  // zero mean, unit-ish variance
+  }
+
+ private:
+  float Uniform() {
+    state_ = state_ * 6364136223846793005ull + 1442695040888963407ull;
+    return static_cast<float>((state_ >> 40) & 0xFFFFFF) / 16777216.0f;
+  }
+  uint64_t state_;
+};
+
+// Fills `mat` with N(0, stddev), and additionally gives a few columns of each
+// row a 10x larger magnitude. Outlier channels are the known hard case for
+// per-tensor int8; per-row/per-column scales are supposed to absorb them.
+void FillNormal(MatStorageT<float>& mat, uint64_t seed, float stddev,
+                bool outliers) {
+  Rng rng(seed);
+  for (size_t r = 0; r < mat.Rows(); ++r) {
+    float* HWY_RESTRICT row = mat.Row(r);
+    for (size_t c = 0; c < mat.Cols(); ++c) {
+      row[c] = rng.Normal() * stddev;
+    }
+    if (outliers) {
+      for (size_t c = (r * 7) % 64; c < mat.Cols(); c += 512) {
+        row[c] *= 10.0f;
+      }
+    }
+    for (size_t c = mat.Cols(); c < mat.Stride(); ++c) row[c] = 0.0f;
+  }
+}
+
+// Converts F32 `in` to `MatT` (BF16 or a compressed stream), row by row.
+template <typename MatT>
+void ConvertRows(const MatStorageT<float>& in, MatStorageT<MatT>& out,
+                 ThreadingContext& ctx) {
+  CompressWorkingSet ws;
+  ws.tls.resize(ctx.pools.MaxWorkers());
+  const size_t cols = in.Cols();
+  ParallelFor(Parallelism::kFlat, in.Rows(), ctx, /*cluster_idx=*/0,
+              Callers::kTest, [&](size_t r, size_t thread) HWY_ATTR {
+                Compress(in.Row(r), cols, ws.tls[thread],
+                         MakeSpan(out.Row(r), cols), /*packed_ofs=*/0);
+              });
+}
+
+//------------------------------------------------------------------------------
+// Reference and error metric
+
+// `B` is transposed: `ref[m, n] = sum_k A[m, k] * B[n, k]`.
+void ReferenceMatMul(const MatStorageT<float>& A, const MatStorageT<float>& B,
+                     MatStorageT<double>& ref, ThreadingContext& ctx) {
+  const size_t K = A.Cols();
+  ParallelFor(Parallelism::kFlat, A.Rows(), ctx, /*cluster_idx=*/0,
+              Callers::kTest, [&](size_t m, size_t /*thread*/) {
+                const float* HWY_RESTRICT a = A.Row(m);
+                double* HWY_RESTRICT out = ref.Row(m);
+                for (size_t n = 0; n < B.Rows(); ++n) {
+                  const float* HWY_RESTRICT b = B.Row(n);
+                  double sum = 0.0;
+                  for (size_t k = 0; k < K; ++k) {
+                    sum += static_cast<double>(a[k]) * static_cast<double>(b[k]);
+                  }
+                  out[n] = sum;
+                }
+              });
+}
+
+// Relative Frobenius error ||C - ref|| / ||ref||.
+template <typename TC>
+double RelError(const MatStorageT<TC>& C, const MatStorageT<double>& ref) {
+  double num = 0.0, den = 0.0;
+  for (size_t r = 0; r < ref.Rows(); ++r) {
+    const TC* HWY_RESTRICT c = C.Row(r);
+    const double* HWY_RESTRICT e = ref.Row(r);
+    for (size_t n = 0; n < ref.Cols(); ++n) {
+      const double d = hwy::ConvertScalarTo<double>(c[n]) - e[n];
+      num += d * d;
+      den += e[n] * e[n];
+    }
+  }
+  return (den == 0.0) ? 0.0 : std::sqrt(num / den);
+}
+
+//------------------------------------------------------------------------------
+// Timing
+
+struct Result {
+  double median_sec = 0.0;
+  double gflops = 0.0;
+  double rel_error = -1.0;  // < 0 if not measured
+};
+
+// Repeats `fn` (which returns the autotune state) until autotuning has settled,
+// then collects `num_samples` timings and returns the median.
+template <class Fn>
+Result TimeMatMul(size_t M, size_t K, size_t N, Fn&& fn) {
+  const size_t num_samples = M < 32 ? 40 : 12;
+  std::vector<double> times;
+  times.reserve(num_samples);
+
+  // Bound the loop: a config that never reports Best() would otherwise hang.
+  // Skip a few runs after autotuning settles, so the first timed sample is not
+  // the one that still has the autotuner's working set in cache.
+  size_t warmup = 3;
+  for (size_t iter = 0; times.size() < num_samples && iter < 8192; ++iter) {
+    const double t0 = hwy::platform::Now();
+    MMPerKey* per_key = fn();
+    const double t1 = hwy::platform::Now();
+    if (!per_key->autotune.Best()) continue;
+    if (warmup != 0) {
+      --warmup;
+      continue;
+    }
+    times.push_back(t1 - t0);
+  }
+  HWY_ASSERT(!times.empty());
+
+  std::sort(times.begin(), times.end());
+  Result r;
+  r.median_sec = times[times.size() / 2];
+  r.gflops = 2.0 * M * K * N / r.median_sec * 1E-9;
+  return r;
+}
+
+//------------------------------------------------------------------------------
+// One shape
+
+// Runs BF16xBF16, BF16xSFP and int8 W8A8 on the same `M x K x N` problem.
+// `check_error` also computes the F64 reference, which is O(M*K*N) scalar work
+// and thus only affordable for smaller shapes.
+void BenchShape(size_t M, size_t K, size_t N, bool check_error,
+                ThreadingContext& ctx, MatMulEnv& env_bf, MatMulEnv& env_sfp,
+                MatMulEnv& env_i8, MMI8AStorage& a_i8,
+                bool a_outliers = true) {
+  const Allocator& allocator = ctx.allocator;
+  const Extents2D A_extents(M, K);
+  const Extents2D B_extents(N, K);  // already transposed
+  const Extents2D C_extents(M, N);
+
+  // Sources, in F32.
+  MatStorageT<float> A_f32("A_f32", A_extents, allocator, MatPadding::kOdd);
+  MatStorageT<float> B_f32("B_f32", B_extents, allocator, MatPadding::kOdd);
+  FillNormal(A_f32, /*seed=*/1, /*stddev=*/1.0f, a_outliers);
+  FillNormal(B_f32, /*seed=*/2, /*stddev=*/0.02f, /*outliers=*/false);
+
+  // Operands for the BF16 kernel.
+  MatStorageT<BF16> A_bf("A_bf", A_extents, allocator, MatPadding::kOdd);
+  MatStorageT<BF16> B_bf("B_bf", B_extents, allocator, MatPadding::kOdd);
+  MatStorageT<SfpStream> B_sfp("B_sfp", B_extents, allocator, MatPadding::kOdd);
+  ConvertRows(A_f32, A_bf, ctx);
+  ConvertRows(B_f32, B_bf, ctx);
+  ConvertRows(B_f32, B_sfp, ctx);
+
+  // Operands for the int8 kernel. `A` is quantized inside `MatMulI8`.
+  MatStorageT<int8_t> B_i8("B_i8", B_extents, allocator, MatPadding::kOdd);
+  hwy::AlignedVector<float> b_scale(N);
+  const MMI8B B_packed = PackB(B_f32, B_i8, b_scale.data(), ctx);
+
+  MatStorageT<float> C_bf("C_bf", C_extents, allocator, MatPadding::kOdd);
+  MatStorageT<float> C_sfp("C_sfp", C_extents, allocator, MatPadding::kOdd);
+  MatStorageT<float> C_i8("C_i8", C_extents, allocator, MatPadding::kOdd);
+  C_bf.AllocateAndAttachRowPtrs(env_bf.row_ptrs);
+  C_sfp.AllocateAndAttachRowPtrs(env_sfp.row_ptrs);
+  C_i8.AllocateAndAttachRowPtrs(env_i8.row_ptrs);
+
+  Tristate use_spinning = Tristate::kDefault;
+  ctx.pools.MaybeStartSpinning(use_spinning);
+
+  const Result r_bf = TimeMatMul(M, K, N, [&] {
+    return MatMul(A_bf, B_bf, /*add=*/nullptr, env_bf, C_bf);
+  });
+  const Result r_sfp = TimeMatMul(M, K, N, [&] {
+    return MatMul(A_bf, B_sfp, /*add=*/nullptr, env_sfp, C_sfp);
+  });
+  const Result r_i8 = TimeMatMul(M, K, N, [&] {
+    return MatMulI8(A_bf, B_packed, /*add=*/nullptr, env_i8, C_i8, a_i8);
+  });
+
+  ctx.pools.MaybeStopSpinning(use_spinning);
+
+  double e_bf = -1.0, e_sfp = -1.0, e_i8 = -1.0;
+  if (check_error) {
+    MatStorageT<double> ref("ref", C_extents, allocator, MatPadding::kOdd);
+    ReferenceMatMul(A_f32, B_f32, ref, ctx);
+    e_bf = RelError(C_bf, ref);
+    e_sfp = RelError(C_sfp, ref);
+    e_i8 = RelError(C_i8, ref);
+  }
+
+  printf("%5zu %6zu %7zu | %8.1f %8.1f %8.1f | %6.3f %6.3f %6.3f | %5.2fx %5.2fx",
+         M, K, N, r_bf.gflops, r_sfp.gflops, r_i8.gflops,
+         r_bf.median_sec * 1E3, r_sfp.median_sec * 1E3, r_i8.median_sec * 1E3,
+         r_i8.gflops / r_bf.gflops, r_i8.gflops / r_sfp.gflops);
+  if (check_error) {
+    printf(" | %.2e %.2e %.2e", e_bf, e_sfp, e_i8);
+  }
+  printf("\n");
+  fflush(stdout);
+}
+
+// Measures raw instruction throughput of the two dot products with 16
+// independent accumulator chains, no memory traffic. The matmul speedups below
+// cannot exceed this ratio, and how close they get says how much of the win is
+// compute rather than the halved footprint of `B`.
+void BenchDotThroughput() {
+  const hn::ScalableTag<BF16> dbf;
+  const hn::Repartition<float, decltype(dbf)> df;
+  const hn::ScalableTag<int8_t> di8;
+  const hn::Repartition<int32_t, decltype(di8)> di32;
+  constexpr size_t kChains = 16;
+  constexpr size_t kReps = 2000000;
+
+  const size_t bf16_macs = hn::Lanes(dbf);  // per instruction
+  const size_t i8_macs = hn::Lanes(di8);
+
+  double keep = 0.0;
+  double bf_sec = 0.0, i8_sec = 0.0;
+
+  {
+    const auto a = hn::Set(dbf, hwy::ConvertScalarTo<BF16>(1.0f));
+    hn::Vec<decltype(df)> c[kChains], unused = hn::Zero(df);
+    for (size_t i = 0; i < kChains; ++i) c[i] = hn::Zero(df);
+    const double t0 = hwy::platform::Now();
+    for (size_t r = 0; r < kReps; ++r) {
+      for (size_t i = 0; i < kChains; ++i) {
+        c[i] = hn::ReorderWidenMulAccumulate(df, a, a, c[i], unused);
+      }
+    }
+    bf_sec = hwy::platform::Now() - t0;
+    for (size_t i = 0; i < kChains; ++i) keep += hn::GetLane(c[i]);
+  }
+  {
+    const auto a = hn::Set(di8, int8_t{1});
+    hn::Vec<decltype(di32)> c[kChains];
+    for (size_t i = 0; i < kChains; ++i) c[i] = hn::Zero(di32);
+    const double t0 = hwy::platform::Now();
+    for (size_t r = 0; r < kReps; ++r) {
+      for (size_t i = 0; i < kChains; ++i) {
+        c[i] = hn::SumOfMulQuadAccumulate(di32, a, a, c[i]);
+      }
+    }
+    i8_sec = hwy::platform::Now() - t0;
+    for (size_t i = 0; i < kChains; ++i) keep += hn::GetLane(c[i]);
+  }
+  hwy::PreventElision(keep);
+
+  const double ops = static_cast<double>(kChains) * kReps;
+  const double bf_gmac = ops * bf16_macs / bf_sec * 1E-9;
+  const double i8_gmac = ops * i8_macs / i8_sec * 1E-9;
+  printf(
+      "1-core dot product throughput: bf16 %.1f GMAC/s, int8 %.1f GMAC/s "
+      "(%.2fx)\n",
+      bf_gmac, i8_gmac, i8_gmac / bf_gmac);
+}
+
+void BenchAll() {
+  ThreadingArgs threading_args;
+  ThreadingContext ctx(threading_args);
+  printf("target=%s %s %s\n", hwy::TargetName(HWY_TARGET),
+         ctx.topology.TopologyString(), ctx.pools.PinString());
+  printf("B biased to u8: %d, HWY_NATIVE_DOT_BF16=%d, vector bytes=%zu\n",
+         GEMMA_MM_I8_BIASED_B, HWY_NATIVE_DOT_BF16,
+         hn::Lanes(hn::ScalableTag<uint8_t>()));
+
+  BenchDotThroughput();
+
+  MatMulEnv env_bf(ctx), env_sfp(ctx), env_i8(ctx);
+  // Sized for the largest shape below.
+  MMI8AStorage a_i8(/*max_M=*/512, /*max_K=*/8192, ctx.allocator);
+
+  printf(
+      "\n    M      K       N |  GFLOPS: bf16   sfp     i8 |     ms: bf16    "
+      "sfp     i8 | i8 vs bf16/sfp | rel err: bf16 sfp i8\n");
+
+  // Gemma3-1B decode shapes, as in `ops/bench_matmul.cc`.
+  for (size_t M : {size_t{1}, size_t{4}}) {
+    BenchShape(M, 1152, 1536, /*check_error=*/false, ctx, env_bf, env_sfp,
+               env_i8, a_i8);  // QKV
+    BenchShape(M, 1152, 13824, false, ctx, env_bf, env_sfp, env_i8,
+               a_i8);  // FFN gate+up
+    BenchShape(M, 6912, 1152, false, ctx, env_bf, env_sfp, env_i8,
+               a_i8);  // FFN down
+    BenchShape(M, 1152, 32768, false, ctx, env_bf, env_sfp, env_i8,
+               a_i8);  // logits (N reduced to fit memory)
+  }
+
+  // Prefill / batched shapes.
+  BenchShape(128, 3072, 3072, false, ctx, env_bf, env_sfp, env_i8, a_i8);
+  BenchShape(512, 3072, 3072, false, ctx, env_bf, env_sfp, env_i8, a_i8);
+  BenchShape(128, 1152, 13824, false, ctx, env_bf, env_sfp, env_i8, a_i8);
+
+  // B larger than last-level cache in both formats, so neither kernel can
+  // hide the streaming cost of B.
+  printf("\n(B exceeds LLC in both formats)\n");
+  BenchShape(1, 4096, 32768, false, ctx, env_bf, env_sfp, env_i8, a_i8);
+  BenchShape(8, 4096, 32768, false, ctx, env_bf, env_sfp, env_i8, a_i8);
+  BenchShape(128, 4096, 32768, false, ctx, env_bf, env_sfp, env_i8, a_i8);
+
+  printf(
+      "\nAccuracy vs F64 reference. A has outlier channels (10x), which is the"
+      "\nknown hard case for per-token int8 activations:\n");
+  BenchShape(32, 1152, 512, /*check_error=*/true, ctx, env_bf, env_sfp, env_i8,
+             a_i8, /*a_outliers=*/true);
+  BenchShape(32, 3072, 512, true, ctx, env_bf, env_sfp, env_i8, a_i8, true);
+  BenchShape(32, 6912, 512, true, ctx, env_bf, env_sfp, env_i8, a_i8, true);
+
+  printf("\nSame, but A is plain Gaussian with no outlier channels:\n");
+  BenchShape(32, 1152, 512, true, ctx, env_bf, env_sfp, env_i8, a_i8,
+             /*a_outliers=*/false);
+  BenchShape(32, 3072, 512, true, ctx, env_bf, env_sfp, env_i8, a_i8, false);
+  BenchShape(32, 6912, 512, true, ctx, env_bf, env_sfp, env_i8, a_i8, false);
+}
+
+// NOLINTNEXTLINE(google-readability-namespace-comments)
+}  // namespace HWY_NAMESPACE
+}  // namespace gcpp
+HWY_AFTER_NAMESPACE();
+
+#if HWY_ONCE
+namespace gcpp {
+HWY_EXPORT(BenchAll);
+}  // namespace gcpp
+
+int main(int /*argc*/, char** /*argv*/) {
+  // Best available target only; this is a benchmark, not a test.
+  HWY_DYNAMIC_DISPATCH(gcpp::BenchAll)();
+  return 0;
+}
+#endif  // HWY_ONCE
diff --git a/ops/matmul-inl.h b/ops/matmul-inl.h
index cabe5392..6cb92d22 100644
--- a/ops/matmul-inl.h
+++ b/ops/matmul-inl.h
@@ -17,6 +17,8 @@
 #include <stdint.h>
 #include <stdio.h>
 
+#include <cmath>
+#include <cstdlib>
 #include <vector>
 
 #include "compression/types.h"
@@ -53,6 +55,30 @@ namespace gcpp {
 namespace HWY_NAMESPACE {
 namespace hn = hwy::HWY_NAMESPACE;
 
+// Shared integer dot-product primitive for direct quantized MatMul kernels.
+// `kWeightsFirst` selects the operand order required by the encoding:
+// unsigned weights must precede signed activations, whereas signed W8A8 keeps
+// the activation first. This keeps target-specific dot-product details out of
+// the W8A8 kernel.
+template <bool kWeightsFirst, class DI32, class VA8, class VB8,
+          class VI32 = hn::Vec<DI32>>
+static HWY_INLINE void MMQuantizedDot4Accumulate(
+    DI32 di32, VA8 a, VB8 b0, VB8 b1, VB8 b2, VB8 b3, VI32& c0, VI32& c1,
+    VI32& c2, VI32& c3) {
+  static_assert(kNR == 4);
+  if constexpr (kWeightsFirst) {
+    c0 = hn::SumOfMulQuadAccumulate(di32, b0, a, c0);
+    c1 = hn::SumOfMulQuadAccumulate(di32, b1, a, c1);
+    c2 = hn::SumOfMulQuadAccumulate(di32, b2, a, c2);
+    c3 = hn::SumOfMulQuadAccumulate(di32, b3, a, c3);
+  } else {
+    c0 = hn::SumOfMulQuadAccumulate(di32, a, b0, c0);
+    c1 = hn::SumOfMulQuadAccumulate(di32, a, b1, c1);
+    c2 = hn::SumOfMulQuadAccumulate(di32, a, b2, c2);
+    c3 = hn::SumOfMulQuadAccumulate(di32, a, b3, c3);
+  }
+}
+
 // Like hn::PromoteOddTo, but uses assembly to avoid an extra vector register.
 template <class DF, class DBF = hn::Repartition<BF16, DF>>
 static hn::VFromD<DF> FastPromoteOddTo(DF df, hn::VFromD<DBF> vbf) {
@@ -229,6 +255,20 @@ class MMStoreHorizontalSumsIntoC {
 // Stateless, wraps member functions.
 class MMDecompress {
  public:
+  // Quality-only experiment requested in #1. When enabled, every activation
+  // row is symmetrically quantized to int8 and immediately dequantized to
+  // BF16 before the existing MatMul. This intentionally adds overhead: its
+  // purpose is to isolate activation-quantization error from a new kernel and
+  // weight quantization.
+  static bool ActivationI8RoundtripEnabled() {
+    static const bool enabled = [] {
+      const char* value = std::getenv("GEMMA_MM_I8_ROUNDTRIP_A");
+      return value != nullptr && value[0] != '\0' &&
+             !(value[0] == '0' && value[1] == '\0');
+    }();
+    return enabled;
+  }
+
   // Decompresses `kNR x kc` from `B[row_b, range_kc.begin()]` to row 0,
   // col 0 of `B_view`. Decompressing SFP is relatively cheap on `AVX3_DL`
   // thanks to its large table lookups, and less so on other targets.
@@ -271,7 +311,9 @@ class MMDecompress {
     if constexpr (IsBF16<TA>()) {
       // We can use a view, regardless of columns/padding, because
       // `MMKernel::LoopKC` supports non-vector multiples.
-      return StridedViewBF(A, 0, 0, A.Cols());
+      const StridedViewBF A_view(A, 0, 0, A.Cols());
+      return ActivationI8RoundtripEnabled() ? RoundtripA(A, A_view, env)
+                                            : A_view;
     } else {
       // Always decompress. To reduce code size/compile time, we no longer
       // support a separate F32 kernel; most A are already BF16. We also only
@@ -279,11 +321,52 @@ class MMDecompress {
       HWY_ASSERT(options.cluster_idx == 0);
       const StridedViewBF A_view = env.A_BF.A(A.Extents());
       AutotuneDecompressA(A, A_view, autotune, env, options);
-      return A_view;
+      return ActivationI8RoundtripEnabled() ? RoundtripA(A, A_view, env)
+                                            : A_view;
     }
   }
 
+  // `TwoMatMul` only accepts BF16 A and does not call `MaybeDecompressA`.
+  static HWY_INLINE StridedViewBF MaybeRoundtripA(const MatPtrT<BF16>& A,
+                                                  const MatMulEnv& env) {
+    const StridedViewBF A_view(A, 0, 0, A.Cols());
+    return ActivationI8RoundtripEnabled() ? RoundtripA(A, A_view, env)
+                                          : A_view;
+  }
+
  private:
+  template <typename TA>
+  static HWY_NOINLINE StridedViewBF RoundtripA(const MatPtrT<TA>& A,
+                                               const StridedViewBF source,
+                                               const MatMulEnv& env) {
+    const StridedViewBF dest = env.A_BF.A(A.Extents());
+    constexpr float kI8Max = 127.0f;
+
+    for (size_t row = 0; row < A.Rows(); ++row) {
+      const BF16* HWY_RESTRICT from = source.Row(row);
+      BF16* HWY_RESTRICT to = dest.Row(row);
+
+      float max_abs = 0.0f;
+      for (size_t col = 0; col < A.Cols(); ++col) {
+        max_abs = HWY_MAX(
+            max_abs,
+            std::fabs(hwy::ConvertScalarTo<float>(from[col])));
+      }
+
+      const float scale = max_abs == 0.0f ? 1.0f : max_abs / kI8Max;
+      const float inv_scale = max_abs == 0.0f ? 0.0f : kI8Max / max_abs;
+      for (size_t col = 0; col < A.Cols(); ++col) {
+        const float value = hwy::ConvertScalarTo<float>(from[col]);
+        const int32_t quantized = HWY_MIN(
+            int32_t{127},
+            HWY_MAX(int32_t{-127}, static_cast<int32_t>(std::lroundf(
+                                            value * inv_scale))));
+        to[col] = hwy::ConvertScalarTo<BF16>(quantized * scale);
+      }
+    }
+    return dest;
+  }
+
   // Decompresses all `M x K` from `A` into padded BF16 `A_view`.
   static HWY_NOINLINE void DecompressA(const MatPtrT<float>& A,
                                        const StridedViewBF A_view,
@@ -395,6 +478,9 @@ class MMDecompress {
 // Stateless, wraps member functions. Contains the innermost 2-4 loops.
 class MMKernel {
  public:
+  // Type of the `A` operand, see `MMLoops::Dispatch`.
+  using AView = StridedViewBF;
+
   // Loop over NC/MC/KC, called from the outer loops. The MOMMS B3A2C0 reads
   // `mc x kc` of A, `nc x kc` of B, and updates the `mc x nc` `C_MC_NC`.
   // `CView` is either `RowPtrs<TC>` or `StridedView<TC>`.
@@ -1217,9 +1303,12 @@ class MMLoops {
  public:
   // Called from `MatMul` from two places: either with the next autotune config,
   // or with the best config. `B2` is null unless called from `TwoMatMul`.
-  template <typename TB, typename TC>
-  static HWY_NOINLINE void Dispatch(const StridedViewBF A, const MatPtrT<TB>& B,
-                                    const MatPtrT<TB>* B2, RowPtrs<TC> C,
+  // `Kernel` is `MMKernel` (BF16) or `MMI8Kernel` (int8); it defines the type
+  // of `A` and, with `BT`, how a tile is computed. The loops themselves only
+  // partition the work, hence they are shared between kernels.
+  template <class Kernel, typename BT, typename TC>
+  static HWY_NOINLINE void Dispatch(const typename Kernel::AView A, const BT& B,
+                                    const BT* B2, RowPtrs<TC> C,
                                     const MMArgs& args) {
     GCPP_ZONE(args.env.ctx, args.env.ctx.Worker(args.options.cluster_idx),
               Zones::kMMDispatch);
@@ -1227,7 +1316,7 @@ class MMLoops {
     DispatchParallelism(
         args.options.parallelism, [&](const auto& parallel) HWY_ATTR {
           DispatchOrder(args.order, [&](const auto& order) HWY_ATTR {
-            Loop(order, parallel, A, B, B2, C, args);
+            Loop<Kernel>(order, parallel, A, B, B2, C, args);
           });
         });
   }
@@ -1240,10 +1329,10 @@ class MMLoops {
   }
 
   // Single M and K ranges, parallel N.
-  template <typename TB, typename TC, class Parallel>
+  template <class Kernel, typename BT, typename TC, class Parallel>
   static HWY_INLINE void Loop(MMOrderNT, Parallel parallel,
-                              const StridedViewBF A, const MatPtrT<TB>& B,
-                              const MatPtrT<TB>* B2, RowPtrs<TC> C,
+                              const typename Kernel::AView A, const BT& B,
+                              const BT* B2, RowPtrs<TC> C,
                               const MMArgs& args) {
     const auto zone = args.env.ctx.profiler_zones.Get(Zones::kMMNT);
     HWY_DASSERT(args.ranges_mc.NumTasks() == 1);
@@ -1258,14 +1347,14 @@ class MMLoops {
           MMZone mm_zone;
           mm_zone.MaybeEnter(worker, zone, args.env, &args.autotune);
 
-          MMKernel::B3A2C0(A, B, range_mc, range_kc, range_nc, args, MMSetC(),
+          Kernel::B3A2C0(A, B, range_mc, range_kc, range_nc, args, MMSetC(),
                            C.View(0, range_nc.begin(), range_nc.Num()));
 
           const StridedViewBF C2 = args.env.C_tiles.C(
               Extents2D(range_mc.Num(), range_nc.Num()), worker);
 
           if (B2 != nullptr) {
-            MMKernel::B3A2C0(A, *B2, range_mc, range_kc, range_nc, args,
+            Kernel::B3A2C0(A, *B2, range_mc, range_kc, range_nc, args,
                              MMSetC(), C2);
           }
 
@@ -1276,10 +1365,10 @@ class MMLoops {
   }
 
   // Single M range, parallel N, sequential K. Sets C, then accumulates.
-  template <typename TB, typename TC, class Parallel>
+  template <class Kernel, typename BT, typename TC, class Parallel>
   static HWY_INLINE void Loop(MMOrderNT_K, Parallel parallel,
-                              const StridedViewBF A, const MatPtrT<TB>& B,
-                              const MatPtrT<TB>* B2, RowPtrs<TC> C,
+                              const typename Kernel::AView A, const BT& B,
+                              const BT* B2, RowPtrs<TC> C,
                               const MMArgs& args) {
     const auto zone = args.env.ctx.profiler_zones.Get(Zones::kMMNT_K);
     HWY_DASSERT(args.ranges_mc.NumTasks() == 1);
@@ -1291,7 +1380,7 @@ class MMLoops {
                   [&](const IndexRange& range_nc, size_t worker) HWY_ATTR {
                     MMZone mm_zone;
                     mm_zone.MaybeEnter(worker, zone, args.env, &args.autotune);
-                    MMKernel::ForeachKC(
+                    Kernel::ForeachKC(
                         A, B, range_mc, args.ranges_kc, range_nc, args,
                         C.View(0, range_nc.begin(), range_nc.Num()));
 
@@ -1299,7 +1388,7 @@ class MMLoops {
                         Extents2D(range_mc.Num(), range_nc.Num()), worker);
 
                     if (B2 != nullptr) {
-                      MMKernel::ForeachKC(A, *B2, range_mc, args.ranges_kc,
+                      Kernel::ForeachKC(A, *B2, range_mc, args.ranges_kc,
                                           range_nc, args, C2);
                     }
 
@@ -1312,10 +1401,10 @@ class MMLoops {
 
   // Parallel loops over mc/nc blocks of M/range_n, single K.
   // Fills `mc x nc` sections of C.
-  template <typename TB, typename TC, class Parallel>
+  template <class Kernel, typename BT, typename TC, class Parallel>
   static HWY_INLINE void Loop(MMOrderNT_MT, Parallel parallel,
-                              const StridedViewBF A, const MatPtrT<TB>& B,
-                              const MatPtrT<TB>* B2, RowPtrs<TC> C,
+                              const typename Kernel::AView A, const BT& B,
+                              const BT* B2, RowPtrs<TC> C,
                               const MMArgs& args) {
     const auto zone = args.env.ctx.profiler_zones.Get(Zones::kMMNT_MT);
     HWY_DASSERT(args.ranges_kc.NumTasks() == 1);
@@ -1327,7 +1416,7 @@ class MMLoops {
             size_t worker) HWY_ATTR {
           MMZone mm_zone;
           mm_zone.MaybeEnter(worker, zone, args.env, &args.autotune);
-          MMKernel::B3A2C0(
+          Kernel::B3A2C0(
               A, B, range_mc, range_kc, range_nc, args, MMSetC(),
               C.View(range_mc.begin(), range_nc.begin(), range_nc.Num()));
 
@@ -1335,7 +1424,7 @@ class MMLoops {
               Extents2D(range_mc.Num(), range_nc.Num()), worker);
 
           if (B2 != nullptr) {
-            MMKernel::B3A2C0(A, *B2, range_mc, range_kc, range_nc, args,
+            Kernel::B3A2C0(A, *B2, range_mc, range_kc, range_nc, args,
                              MMSetC(), C2);
           }
           if constexpr (IsBF16<TC>()) {
@@ -1346,10 +1435,10 @@ class MMLoops {
 
   // Parallel loops over mc/nc blocks of M/range_n, sequential K.
   // Accumulates into `mc x nc` sections of `C`.
-  template <typename TB, typename TC, class Parallel>
+  template <class Kernel, typename BT, typename TC, class Parallel>
   static HWY_INLINE void Loop(MMOrderNT_MT_K, Parallel parallel,
-                              const StridedViewBF A, const MatPtrT<TB>& B,
-                              const MatPtrT<TB>* B2, RowPtrs<TC> C,
+                              const typename Kernel::AView A, const BT& B,
+                              const BT* B2, RowPtrs<TC> C,
                               const MMArgs& args) {
     const auto zone = args.env.ctx.profiler_zones.Get(Zones::kMMNT_MT_K);
 
@@ -1359,7 +1448,7 @@ class MMLoops {
             size_t worker) HWY_ATTR {
           MMZone mm_zone;
           mm_zone.MaybeEnter(worker, zone, args.env, &args.autotune);
-          MMKernel::ForeachKC(
+          Kernel::ForeachKC(
               A, B, range_mc, args.ranges_kc, range_nc, args,
               C.View(range_mc.begin(), range_nc.begin(), range_nc.Num()));
 
@@ -1367,7 +1456,7 @@ class MMLoops {
               Extents2D(range_mc.Num(), range_nc.Num()), worker);
 
           if (B2 != nullptr) {
-            MMKernel::ForeachKC(A, *B2, range_mc, args.ranges_kc, range_nc,
+            Kernel::ForeachKC(A, *B2, range_mc, args.ranges_kc, range_nc,
                                 args, C2);
           }
 
@@ -1378,10 +1467,10 @@ class MMLoops {
   }
 
   // Parallel loops over mc/nc blocks of M/range_n via SFC, single K.
-  template <typename TB, typename TC, class Parallel>
+  template <class Kernel, typename BT, typename TC, class Parallel>
   static HWY_INLINE void Loop(MMOrderSFC, Parallel parallel,
-                              const StridedViewBF A, const MatPtrT<TB>& B,
-                              const MatPtrT<TB>* B2, RowPtrs<TC> C,
+                              const typename Kernel::AView A, const BT& B,
+                              const BT* B2, RowPtrs<TC> C,
                               const MMArgs& args) {
     const auto zone = args.env.ctx.profiler_zones.Get(Zones::kMMSFC);
     HWY_DASSERT(args.ranges_kc.NumTasks() == 1);
@@ -1393,7 +1482,7 @@ class MMLoops {
             size_t worker) HWY_ATTR {
           MMZone mm_zone;
           mm_zone.MaybeEnter(worker, zone, args.env, &args.autotune);
-          MMKernel::B3A2C0(
+          Kernel::B3A2C0(
               A, B, range_mc, range_kc, range_nc, args, MMSetC(),
               C.View(range_mc.begin(), range_nc.begin(), range_nc.Num()));
 
@@ -1401,7 +1490,7 @@ class MMLoops {
               Extents2D(range_mc.Num(), range_nc.Num()), worker);
 
           if (B2 != nullptr) {
-            MMKernel::B3A2C0(A, *B2, range_mc, range_kc, range_nc, args,
+            Kernel::B3A2C0(A, *B2, range_mc, range_kc, range_nc, args,
                              MMSetC(), C2);
           }
           if constexpr (IsBF16<TC>()) {
@@ -1411,10 +1500,10 @@ class MMLoops {
   }
 
   // Parallel loops over mc/nc blocks of M/range_n via SFC, sequential K.
-  template <typename TB, typename TC, class Parallel>
+  template <class Kernel, typename BT, typename TC, class Parallel>
   static HWY_INLINE void Loop(MMOrderSFC_K, Parallel parallel,
-                              const StridedViewBF A, const MatPtrT<TB>& B,
-                              const MatPtrT<TB>* B2, RowPtrs<TC> C,
+                              const typename Kernel::AView A, const BT& B,
+                              const BT* B2, RowPtrs<TC> C,
                               const MMArgs& args) {
     const auto zone = args.env.ctx.profiler_zones.Get(Zones::kMMSFC_K);
 
@@ -1424,7 +1513,7 @@ class MMLoops {
             size_t worker) HWY_ATTR {
           MMZone mm_zone;
           mm_zone.MaybeEnter(worker, zone, args.env, &args.autotune);
-          MMKernel::ForeachKC(
+          Kernel::ForeachKC(
               A, B, range_mc, args.ranges_kc, range_nc, args,
               C.View(range_mc.begin(), range_nc.begin(), range_nc.Num()));
 
@@ -1432,7 +1521,7 @@ class MMLoops {
               Extents2D(range_mc.Num(), range_nc.Num()), worker);
 
           if (B2 != nullptr) {
-            MMKernel::ForeachKC(A, *B2, range_mc, args.ranges_kc, range_nc,
+            Kernel::ForeachKC(A, *B2, range_mc, args.ranges_kc, range_nc,
                                 args, C2);
           }
 
@@ -1486,7 +1575,8 @@ HWY_NOINLINE MMPerKey* MatMul(const MatPtrT<TA>& A, const MatPtrT<TB>& B,
   // BRGeMM path for BF16×BF16 on Intel AMX/AVX-512.
   // Requires M,N,K >= 32 and K % 32 == 0 (AMX tile constraint).
   if constexpr (IsBF16<TA>() && IsBF16<TB>()) {
-    if (M >= 32 && N >= 32 && K >= 32 && (K % 32) == 0) {
+    if (!MMDecompress::ActivationI8RoundtripEnabled() && M >= 32 && N >= 32 &&
+        K >= 32 && (K % 32) == 0) {
       const float scale = A.Scale() * B.Scale();
       MMAutoTune<BRGeMMConfig>& brg_tuner = per_key.brgemm_autotune;
 
@@ -1528,7 +1618,7 @@ HWY_NOINLINE MMPerKey* MatMul(const MatPtrT<TA>& A, const MatPtrT<TB>& B,
   // OneDNN matmul-primitive path for BF16xBF16 via the threadpool runtime.
   // M == 1 was showing worse performance with OneDNN.
   if constexpr (IsBF16<TA>() && IsBF16<TB>()) {
-    if (M > 1) {
+    if (!MMDecompress::ActivationI8RoundtripEnabled() && M > 1) {
       const float scale = A.Scale() * B.Scale();
       if (DoMatMul_OneDnn(A, B, C_rows, M, K, N, scale, add, env,
                           cluster_idx)) {
@@ -1550,7 +1640,7 @@ HWY_NOINLINE MMPerKey* MatMul(const MatPtrT<TA>& A, const MatPtrT<TB>& B,
   if (HWY_LIKELY(tuner.Best())) {
     const MMArgs args(env, M, K, N, A.Scale(), add, options, tuner,
                       *tuner.Best());
-    MMLoops::Dispatch(A_view, B, B2, C_rows, args);
+    MMLoops::Dispatch<MMKernel>(A_view, B, B2, C_rows, args);
     return &per_key;
   }
 
@@ -1570,7 +1660,7 @@ HWY_NOINLINE MMPerKey* MatMul(const MatPtrT<TA>& A, const MatPtrT<TB>& B,
   const MMArgs args(env, M, K, N, A.Scale(), add, options, tuner, cfg);
 
   const uint64_t t0 = hwy::timer::Start();
-  MMLoops::Dispatch(A_view, B, B2, C_rows, args);
+  MMLoops::Dispatch<MMKernel>(A_view, B, B2, C_rows, args);
   MMImpl::NotifyAutotuneResult(env, M, K, N, num_B, t0, tuner, cfg);
 
   return &per_key;
@@ -1603,14 +1693,14 @@ HWY_NOINLINE MMPerKey* TwoMatMul(const MatPtrT<BF16>& A, const MatPtrT<TB>& B1,
       M, K, N, num_B, cache.VectorBytes(), env.per_cluster[cluster_idx]);
 
   // (Also auto-tunes, hence outside the timed section to prevent interference.)
-  const StridedViewBF A_view(A, 0, 0, A.Cols());
+  const StridedViewBF A_view = MMDecompress::MaybeRoundtripA(A, env);
 
   MMAutoTune<MMConfig>& tuner = per_key.autotune;
   if (HWY_LIKELY(tuner.Best())) {
     // Only A scale - B1/B2 may differ, and are passed separately.
     const MMArgs args(env, M, K, N, A.Scale(),
                       /*add=*/nullptr, options, tuner, *tuner.Best());
-    MMLoops::Dispatch(A_view, B1, &B2, C_rows, args);
+    MMLoops::Dispatch<MMKernel>(A_view, B1, &B2, C_rows, args);
     return &per_key;
   }
 
@@ -1634,7 +1724,7 @@ HWY_NOINLINE MMPerKey* TwoMatMul(const MatPtrT<BF16>& A, const MatPtrT<TB>& B1,
                     cfg);
 
   const uint64_t t0 = hwy::timer::Start();
-  MMLoops::Dispatch(A_view, B1, &B2, C_rows, args);
+  MMLoops::Dispatch<MMKernel>(A_view, B1, &B2, C_rows, args);
   MMImpl::NotifyAutotuneResult(env, M, K, N, num_B, t0, tuner, cfg);
 
   return &per_key;
diff --git a/ops/matmul_i8-inl.h b/ops/matmul_i8-inl.h
new file mode 100644
index 00000000..911fb0ef
--- /dev/null
+++ b/ops/matmul_i8-inl.h
@@ -0,0 +1,862 @@
+// Copyright 2025 Google LLC
+// SPDX-License-Identifier: Apache-2.0
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// W8A8 MatMul: symmetric int8 weights times symmetric int8 activations,
+// accumulating in int32 via the 4-way dot product (`vpdpbusd` on x86 VNNI,
+// `sdot`/`usdot` on NEON, `svdot` on SVE). Unlike `MatMul`, which dequantizes
+// `B` to BF16 for every tile (see `MMDecompress::DecompressB`), `B` is
+// consumed as-is.
+//
+// Quantization scheme (see `#560` discussion): per-row (per-token) scales for
+// `A`, computed on the fly, and per-row-of-transposed-B (per output channel)
+// scales baked in at pack time. Both are symmetric, i.e. no zero point, so
+// `C[r, c] = a_scale[r] * b_scale[c] * dot(qa[r], qb[c])` and the int32
+// accumulation can run over an entire `kc` range before a single scaling step.
+//
+// On x86 the 4-way dot product requires one unsigned operand, so `B` is biased
+// by 128 and the `128 * sum_k(qa)` term is subtracted per `kc` range, using
+// prefix sums of the quantized `A`. Biasing `B` rather than `A` is what makes
+// that per-range correction cheap: the correction then depends on `A`, which is
+// small and quantized per call anyway, instead of on `B`. It also keeps the
+// values written to `C` close to the true partial sums; correcting once over
+// the whole `K` would inflate the intermediates that `MMAddC` accumulates
+// through `C`, which loses a lot of precision when `C` is BF16 and the weight
+// channels are not zero-mean.
+
+#include <stddef.h>
+#include <stdint.h>
+
+#include <cmath>
+#include <cstdlib>
+
+#include "ops/matmul.h"  // IWYU pragma: export
+#include "util/basics.h"
+#include "util/mat.h"
+#include "hwy/base.h"
+
+// Include guard for (potentially) SIMD code.
+#if defined(THIRD_PARTY_GEMMA_CPP_MATMUL_I8_TOGGLE) == \
+    defined(HWY_TARGET_TOGGLE)
+#ifdef THIRD_PARTY_GEMMA_CPP_MATMUL_I8_TOGGLE
+#undef THIRD_PARTY_GEMMA_CPP_MATMUL_I8_TOGGLE
+#else
+#define THIRD_PARTY_GEMMA_CPP_MATMUL_I8_TOGGLE
+#endif
+
+#include "hwy/highway.h"
+// After highway.h
+#include "compression/compress-inl.h"
+#include "ops/matmul-inl.h"
+
+// `SumOfMulQuadAccumulate` is native for i8*i8 on NEON with `FEAT_DotProd`
+// and on SVE, but on x86 only for u8*i8 (`vpdpbusd`); there, i8*i8 costs two
+// VNNI ops plus a shift and subtract, which would give up most of the win.
+// Hence bias `B` by 128 into u8 on x86, and correct for it via `A`.
+// Define `GEMMA_MM_I8_FORCE_BIASED_B` to 0 or 1 to exercise either encoding
+// regardless of target; both are correct everywhere, only the speed differs.
+// `ops/matmul_i8_test.cc` is built twice, once each way.
+#undef GEMMA_MM_I8_BIASED_B
+#ifdef GEMMA_MM_I8_FORCE_BIASED_B
+#define GEMMA_MM_I8_BIASED_B GEMMA_MM_I8_FORCE_BIASED_B
+#elif HWY_TARGET <= HWY_AVX3_DL
+#define GEMMA_MM_I8_BIASED_B 1
+#else
+#define GEMMA_MM_I8_BIASED_B 0
+#endif
+
+HWY_BEFORE_NAMESPACE();
+namespace gcpp {
+namespace HWY_NAMESPACE {
+namespace hn = hwy::HWY_NAMESPACE;
+
+// `A` is always symmetric int8; only `B`'s encoding varies by target, see
+// `GEMMA_MM_I8_BIASED_B`.
+using MMI8AT = int8_t;
+#if GEMMA_MM_I8_BIASED_B
+using MMI8BT = uint8_t;
+#else
+using MMI8BT = int8_t;
+#endif
+
+// Largest quantized magnitude. 127 rather than 128 keeps the scheme symmetric,
+// which is what lets us skip the zero-point correction terms.
+HWY_INLINE_VAR constexpr float kMMI8Max = 127.0f;
+
+// Experimental QuaRot-style preprocessing. Applying the same orthonormal
+// transform to A and each row of transposed B leaves their dot product
+// unchanged, while spreading isolated activation outliers over a block. The
+// fixed signs avoid always applying the same Hadamard basis to every block.
+// 128 divides all Gemma 3 1B MatMul K dimensions.
+HWY_INLINE_VAR constexpr size_t kMMI8RotateBlock = 128;
+
+static inline bool MMI8RotateEnabled() {
+  static const bool enabled = [] {
+    const char* value = std::getenv("GEMMA_MM_I8_ROTATE");
+    return value != nullptr && value[0] != '\0' &&
+           !(value[0] == '0' && value[1] == '\0');
+  }();
+  return enabled;
+}
+
+static inline bool MMI8CanRotate(size_t k) {
+  return MMI8RotateEnabled() && (k % kMMI8RotateBlock) == 0;
+}
+
+static HWY_NOINLINE void MMI8Rotate(float* HWY_RESTRICT row, size_t k) {
+  HWY_DASSERT((k % kMMI8RotateBlock) == 0);
+  constexpr float kNormalize = 0.08838834764831845f;  // 1 / sqrt(128)
+  for (size_t block = 0; block < k; block += kMMI8RotateBlock) {
+    float* HWY_RESTRICT x = row + block;
+    for (size_t i = 0; i < kMMI8RotateBlock; ++i) {
+      // Deterministic Rademacher diagonal, shared by A and B.
+      const uint32_t hash =
+          static_cast<uint32_t>(block + i) * 0x9E3779B9u + 0x7F4A7C15u;
+      if ((hash >> 31) != 0) x[i] = -x[i];
+    }
+    for (size_t width = 1; width < kMMI8RotateBlock; width *= 2) {
+      for (size_t start = 0; start < kMMI8RotateBlock; start += 2 * width) {
+        for (size_t i = 0; i < width; ++i) {
+          const float left = x[start + i];
+          const float right = x[start + width + i];
+          x[start + i] = left + right;
+          x[start + width + i] = left - right;
+        }
+      }
+    }
+    for (size_t i = 0; i < kMMI8RotateBlock; ++i) x[i] *= kNormalize;
+  }
+}
+
+//------------------------------------------------------------------------------
+// Quantized operands
+
+// View into quantized `A`, analogous to `StridedViewBF` but carrying the
+// per-row scales (and, when `B` is biased, prefix sums along `K`) alongside,
+// because `MMLoops` passes only this one object down to the kernel.
+struct MMI8AView {
+  // Returns 2D subrange whose top-left is `r, c`, as `StridedView::View`.
+  // Only called on the whole-matrix view, hence the offsets do not compound.
+  MMI8AView View(size_t r, size_t c, size_t cols) const {
+    return MMI8AView{data.View(r, c, cols), scale + r,
+                     prefix + r * prefix_stride + c, prefix_stride};
+  }
+
+  // Sum of the quantized values of row `r` over the `cols` columns of this
+  // view. `prefix` has `K + 1` entries per row, so this is exact for any range.
+  int32_t RowSum(size_t r, size_t cols) const {
+    const int32_t* HWY_RESTRICT p = prefix + r * prefix_stride;
+    return p[cols] - p[0];
+  }
+
+  StridedView<MMI8AT> data;
+  const float* HWY_RESTRICT scale;   // one per row of `data`
+  const int32_t* HWY_RESTRICT prefix;  // null unless `GEMMA_MM_I8_BIASED_B`
+  size_t prefix_stride;
+};
+
+// Transposed, symmetric-int8 `B`: `N` rows of `K` values each, so that a
+// row of `B` is contiguous along `K` and thus already in the layout the 4-way
+// dot product wants. The stored bytes are `q + 128` if `GEMMA_MM_I8_BIASED_B`,
+// else `q`; the buffer is typed `int8_t` either way and reinterpreted in the
+// kernel. Production would pick one encoding for the on-disk format rather
+// than deriving it from the target.
+struct MMI8B {
+  size_t Rows() const { return data->Rows(); }
+  size_t Cols() const { return data->Cols(); }
+
+  const MatPtrT<int8_t>* data;
+  const float* HWY_RESTRICT scale;  // [N] dequantization scale
+};
+
+//------------------------------------------------------------------------------
+// Reduction and store
+
+// Like `MMStoreHorizontalSumsIntoC`, but the tile accumulators are int32 and
+// the scale is a per-row times per-column outer product rather than a scalar.
+template <size_t kRowsAC>
+class MMI8StoreHorizontalSumsIntoC {
+ public:
+  static_assert(kNR == 4);  // for `StoreInterleaved4`
+
+  // Horizontal sums of the 16 (`kRowsAC x kNR`) int32 accumulators, using the
+  // same vector-length-agnostic transpose as the BF16 kernel. Valid because
+  // the 4-way dot product, like BF16's pairwise add, only permutes the terms
+  // of each dot product and thus preserves the horizontal sum.
+  template <class DI32, class VI32 = hn::Vec<DI32>,
+            class D4 = hn::Full128<int32_t>, class V4 = hn::Vec<D4>>
+  HWY_INLINE void Reduce4x4(DI32 di32,                              //
+                            VI32 C00, VI32 C01, VI32 C02, VI32 C03,  //
+                            VI32 C10, VI32 C11, VI32 C12, VI32 C13,  //
+                            VI32 C20, VI32 C21, VI32 C22, VI32 C23,  //
+                            VI32 C30, VI32 C31, VI32 C32, VI32 C33,  //
+                            V4& sum0, V4& sum1, V4& sum2, V4& sum3) {
+    HWY_ALIGN int32_t buf[16 * hn::MaxLanes(di32)];
+    HWY_LANES_CONSTEXPR const size_t N = hn::Lanes(di32);
+
+    MaybeStoreInterleaved4<0>(di32, N, C00, C01, C02, C03, buf);
+    MaybeStoreInterleaved4<1>(di32, N, C10, C11, C12, C13, buf);
+    MaybeStoreInterleaved4<2>(di32, N, C20, C21, C22, C23, buf);
+    MaybeStoreInterleaved4<3>(di32, N, C30, C31, C32, C33, buf);
+
+    const D4 d4;
+    sum0 = MaybeLoad<0>(d4, N, buf);
+    sum1 = MaybeLoad<1>(d4, N, buf);
+    sum2 = MaybeLoad<2>(d4, N, buf);
+    sum3 = MaybeLoad<3>(d4, N, buf);
+
+    for (size_t lane = 1; lane < N; ++lane) {
+      sum0 = MaybeAdd<0>(d4, N, sum0, buf + kNR * lane);
+      sum1 = MaybeAdd<1>(d4, N, sum1, buf + kNR * lane);
+      sum2 = MaybeAdd<2>(d4, N, sum2, buf + kNR * lane);
+      sum3 = MaybeAdd<3>(d4, N, sum3, buf + kNR * lane);
+    }
+  }
+
+  // Dequantizes the four 4-wide int32 dot products and stores them to `C`.
+  // `b_scale` points to the `kNR` current columns and `a_scale` to the current
+  // `range_mc` (hence indexed by `imc + kRow`), whereas `a_rowsum` holds just
+  // this tile's `kRowsAC` values and is indexed by `kRow` alone. It is the sum
+  // of the quantized `A` values over this `kc` range, which undoes `B`'s 128
+  // bias, and is unused when `B` is not biased.
+  template <class D4I, class V4I = hn::Vec<D4I>, class Tag, class CView>
+  HWY_INLINE void Store(D4I d4i, V4I sum0, V4I sum1, V4I sum2, V4I sum3,
+                        const float* HWY_RESTRICT a_scale,
+                        const int32_t* HWY_RESTRICT a_rowsum,
+                        const float* HWY_RESTRICT b_scale,
+                        const float* HWY_RESTRICT add, const size_t imc,
+                        Tag tag, CView C_MC_NR) const {
+    const hn::Full128<float> d4;
+    using V4F = hn::Vec<decltype(d4)>;
+
+    const V4F vb_scale = hn::LoadU(d4, b_scale);
+    HWY_ALIGN static constexpr float kZero[4] = {};
+    const V4F vadd = hn::Load(d4, add ? add : kZero);
+
+    // Each term is `(qb + 128) * qa` instead of `qb * qa`, hence subtract
+    // `128 * sum_k(qa)` over this `kc` range. Applied on every visit, so the
+    // values written to `C` stay close to the true partial sums.
+    MaybeScaleAndStore<0>(d4i, d4, sum0, a_rowsum, vb_scale, vadd, a_scale, tag,
+                          imc, C_MC_NR);
+    MaybeScaleAndStore<1>(d4i, d4, sum1, a_rowsum, vb_scale, vadd, a_scale, tag,
+                          imc, C_MC_NR);
+    MaybeScaleAndStore<2>(d4i, d4, sum2, a_rowsum, vb_scale, vadd, a_scale, tag,
+                          imc, C_MC_NR);
+    MaybeScaleAndStore<3>(d4i, d4, sum3, a_rowsum, vb_scale, vadd, a_scale, tag,
+                          imc, C_MC_NR);
+  }
+
+ private:
+  template <size_t kRow, class DI32, class VI32 = hn::Vec<DI32>>
+  static HWY_INLINE void MaybeStoreInterleaved4(DI32 di32, size_t N, VI32 Cr0,
+                                                VI32 Cr1, VI32 Cr2, VI32 Cr3,
+                                                int32_t* HWY_RESTRICT buf) {
+    if constexpr (kRow < kRowsAC) {
+      hn::StoreInterleaved4(Cr0, Cr1, Cr2, Cr3, di32, buf + 4 * kRow * N);
+    }
+  }
+
+  template <size_t kRow, class D4I, class V4I = hn::Vec<D4I>>
+  static HWY_INLINE V4I MaybeLoad(D4I d4i, size_t N,
+                                  const int32_t* HWY_RESTRICT buf) {
+    if constexpr (kRow < kRowsAC) {
+      return hn::Load(d4i, buf + 4 * kRow * N);
+    } else {
+      return hn::Zero(d4i);
+    }
+  }
+
+  template <size_t kRow, class D4I, class V4I = hn::Vec<D4I>>
+  static HWY_INLINE V4I MaybeAdd(D4I d4i, size_t N, V4I sum,
+                                 const int32_t* HWY_RESTRICT buf) {
+    if constexpr (kRow < kRowsAC) {
+      return hn::Add(sum, hn::Load(d4i, buf + 4 * kRow * N));
+    } else {
+      return sum;
+    }
+  }
+
+  template <size_t kRow, /*deduced:*/ class D4I, class V4I = hn::Vec<D4I>,
+            class D4F, class V4F = hn::Vec<D4F>, class Tag, class CView>
+  static HWY_INLINE void MaybeScaleAndStore(
+      D4I d4i, D4F d4, V4I sum, const int32_t* HWY_RESTRICT a_rowsum,
+      V4F vb_scale, V4F vadd, const float* HWY_RESTRICT a_scale, Tag,
+      const size_t imc, CView C_MC_NR) {
+    if constexpr (kRow < kRowsAC) {
+      using TC = hwy::RemoveCvRef<decltype(C_MC_NR.Row(0)[0])>;
+      TC* HWY_RESTRICT pos = C_MC_NR.Row(imc + kRow);
+      const hn::Rebind<TC, D4F> dc4;
+
+      const V4F vscale = hn::Mul(vb_scale, hn::Set(d4, a_scale[imc + kRow]));
+      if constexpr (GEMMA_MM_I8_BIASED_B) {
+        sum = hn::Sub(sum, hn::Set(d4i, static_cast<int32_t>(
+                                         a_rowsum[kRow] * 128)));
+      }
+      const V4F dot = hn::ConvertTo(d4, sum);
+
+      if constexpr (hwy::IsSame<Tag, MMAddC>()) {
+        vadd = F32FromTC(dc4, hn::Load(dc4, pos));  // load prior value
+      } else {
+        static_assert(hwy::IsSame<Tag, MMSetC>());
+        // vadd remains the bias (added once, the first time we store to C)
+      }
+      const V4F out = hn::MulAdd(dot, vscale, vadd);
+      hn::Store(TCFromF32(dc4, out), dc4, pos);
+    }
+  }
+};  // MMI8StoreHorizontalSumsIntoC
+
+//------------------------------------------------------------------------------
+// Kernel
+
+// Drop-in replacement for `MMKernel` (same `B3A2C0`/`ForeachKC` interface, so
+// that `MMLoops` can drive either), but with int8 operands.
+class MMI8Kernel {
+ public:
+  using AView = MMI8AView;
+
+  template <typename BT, typename Tag, class CView>
+  static void B3A2C0(const AView A, const BT& B, const IndexRange& range_mc,
+                     const IndexRange& range_kc, const IndexRange& range_nc,
+                     const MMArgs& args, Tag out_tag, CView C_MC_NC) {
+    const size_t kc = range_kc.Num();
+    const AView A_view = A.View(range_mc.begin(), range_kc.begin(), kc);
+
+    for (size_t inc = 0; inc < range_nc.Num(); inc += kNR) {
+      // For `add` and `B`, which are global, unlike `C_MC_NC`.
+      const size_t row_b = range_nc.begin() + inc;
+      // No decompression: `B` is already in the layout the kernel wants.
+      const StridedView<int8_t> B_view(*B.data, row_b, range_kc.begin(), kc);
+      const CView C_MC_NR = C_MC_NC.View(0, inc, kNR);
+      const float* HWY_RESTRICT add = args.add ? args.add + row_b : nullptr;
+      A2C0(A_view, B_view, B.scale + row_b, args.mr, range_mc, kc, add, out_tag,
+           C_MC_NR);
+    }
+  }
+
+  template <typename BT, class CView>
+  static void ForeachKC(const AView A, const BT& B, const IndexRange& range_mc,
+                        const IndexRangePartition& ranges_kc,
+                        const IndexRange& range_nc, const MMArgs& args,
+                        CView C_MC_NC) {
+    ranges_kc.VisitFirst([&](const IndexRange& range_kc) {
+      B3A2C0(A, B, range_mc, range_kc, range_nc, args, MMSetC(), C_MC_NC);
+    });
+    ranges_kc.VisitRemaining([&](const IndexRange& range_kc) {
+      B3A2C0(A, B, range_mc, range_kc, range_nc, args, MMAddC(), C_MC_NC);
+    });
+  }
+
+ private:
+  // Innermost loop over `kc` columns in steps of one int8 vector, for
+  // `kRowsAC` rows of `A_view` and `kNR` rows of `B_view`. Mirrors
+  // `MMKernel::LoopKC`: elementwise along `K` with 16 accumulators whose
+  // horizontal sums are the `kRowsAC x kNR` results.
+  template <size_t kRowsAC, /*deduced:*/ class Tag, class CView>
+  static HWY_INLINE void LoopKC(const AView A_view,
+                                const StridedView<int8_t> B_view,
+                                const float* HWY_RESTRICT b_scale, size_t imc,
+                                size_t kc, const float* HWY_RESTRICT add,
+                                Tag tag, CView C_MC_NR) {
+    const hn::ScalableTag<MMI8AT> da8;  // A: always i8
+    const hn::ScalableTag<MMI8BT> db8;  // B: u8 or i8, same lane count
+    const hn::Repartition<int32_t, decltype(da8)> di32;
+    using VA8 = hn::Vec<decltype(da8)>;
+    using VB8 = hn::Vec<decltype(db8)>;
+    using VI32 = hn::Vec<decltype(di32)>;
+    HWY_LANES_CONSTEXPR const size_t N8 = hn::Lanes(da8);
+
+    HWY_DASSERT(kRowsAC <= kMaxMR);
+    static_assert(kNR == 4);
+
+    const MMI8AT* HWY_RESTRICT ar0 = A_view.data.Row(imc + 0);
+    const MMI8AT* HWY_RESTRICT ar1 =
+        kRowsAC > 1 ? A_view.data.Row(imc + 1) : nullptr;
+    const MMI8AT* HWY_RESTRICT ar2 =
+        kRowsAC > 2 ? A_view.data.Row(imc + 2) : nullptr;
+    const MMI8AT* HWY_RESTRICT ar3 =
+        kRowsAC > 3 ? A_view.data.Row(imc + 3) : nullptr;
+    const MMI8BT* HWY_RESTRICT br0 =
+        HWY_RCAST_ALIGNED(const MMI8BT*, B_view.Row(0));
+    const MMI8BT* HWY_RESTRICT br1 =
+        HWY_RCAST_ALIGNED(const MMI8BT*, B_view.Row(1));
+    const MMI8BT* HWY_RESTRICT br2 =
+        HWY_RCAST_ALIGNED(const MMI8BT*, B_view.Row(2));
+    const MMI8BT* HWY_RESTRICT br3 =
+        HWY_RCAST_ALIGNED(const MMI8BT*, B_view.Row(3));
+
+    VI32 C00 = hn::Zero(di32), C01 = hn::Zero(di32), C02 = hn::Zero(di32),
+         C03 = hn::Zero(di32), C10 = hn::Zero(di32), C11 = hn::Zero(di32),
+         C12 = hn::Zero(di32), C13 = hn::Zero(di32), C20 = hn::Zero(di32),
+         C21 = hn::Zero(di32), C22 = hn::Zero(di32), C23 = hn::Zero(di32),
+         C30 = hn::Zero(di32), C31 = hn::Zero(di32), C32 = hn::Zero(di32),
+         C33 = hn::Zero(di32);
+
+    size_t ikc = 0;
+    if (kc >= N8) {
+      HWY_UNROLL(1)
+      for (; ikc <= kc - N8; ikc += N8) {
+        const VB8 b0 = hn::LoadU(db8, br0 + ikc);
+        const VB8 b1 = hn::LoadU(db8, br1 + ikc);
+        const VB8 b2 = hn::LoadU(db8, br2 + ikc);
+        const VB8 b3 = hn::LoadU(db8, br3 + ikc);
+
+        {
+          const VA8 a0 = hn::LoadU(da8, ar0 + ikc);
+          MMQuantizedDot4Accumulate<GEMMA_MM_I8_BIASED_B>(
+              di32, a0, b0, b1, b2, b3, C00, C01, C02, C03);
+        }
+        if constexpr (kRowsAC > 1) {
+          const VA8 a1 = hn::LoadU(da8, ar1 + ikc);
+          MMQuantizedDot4Accumulate<GEMMA_MM_I8_BIASED_B>(
+              di32, a1, b0, b1, b2, b3, C10, C11, C12, C13);
+        }
+        if constexpr (kRowsAC > 2) {
+          const VA8 a2 = hn::LoadU(da8, ar2 + ikc);
+          MMQuantizedDot4Accumulate<GEMMA_MM_I8_BIASED_B>(
+              di32, a2, b0, b1, b2, b3, C20, C21, C22, C23);
+        }
+        if constexpr (kRowsAC > 3) {
+          const VA8 a3 = hn::LoadU(da8, ar3 + ikc);
+          MMQuantizedDot4Accumulate<GEMMA_MM_I8_BIASED_B>(
+              di32, a3, b0, b1, b2, b3, C30, C31, C32, C33);
+        }
+      }
+    }
+
+    // Remainder. `LoadN` zeroes the upper lanes of both operands, so their
+    // products are zero. Zeroing `A` is what makes this safe: a zero `B` lane
+    // does not mean zero in the biased-u8 encoding.
+    const size_t remaining_kc = kc - ikc;
+    HWY_DASSERT(remaining_kc < N8);
+    if (HWY_UNLIKELY(remaining_kc != 0)) {
+      const VB8 b0 = hn::LoadN(db8, br0 + ikc, remaining_kc);
+      const VB8 b1 = hn::LoadN(db8, br1 + ikc, remaining_kc);
+      const VB8 b2 = hn::LoadN(db8, br2 + ikc, remaining_kc);
+      const VB8 b3 = hn::LoadN(db8, br3 + ikc, remaining_kc);
+
+      {
+        const VA8 a0 = hn::LoadN(da8, ar0 + ikc, remaining_kc);
+        MMQuantizedDot4Accumulate<GEMMA_MM_I8_BIASED_B>(
+            di32, a0, b0, b1, b2, b3, C00, C01, C02, C03);
+      }
+      if constexpr (kRowsAC > 1) {
+        const VA8 a1 = hn::LoadN(da8, ar1 + ikc, remaining_kc);
+        MMQuantizedDot4Accumulate<GEMMA_MM_I8_BIASED_B>(
+            di32, a1, b0, b1, b2, b3, C10, C11, C12, C13);
+      }
+      if constexpr (kRowsAC > 2) {
+        const VA8 a2 = hn::LoadN(da8, ar2 + ikc, remaining_kc);
+        MMQuantizedDot4Accumulate<GEMMA_MM_I8_BIASED_B>(
+            di32, a2, b0, b1, b2, b3, C20, C21, C22, C23);
+      }
+      if constexpr (kRowsAC > 3) {
+        const VA8 a3 = hn::LoadN(da8, ar3 + ikc, remaining_kc);
+        MMQuantizedDot4Accumulate<GEMMA_MM_I8_BIASED_B>(
+            di32, a3, b0, b1, b2, b3, C30, C31, C32, C33);
+      }
+    }
+
+    // Sums of the quantized `A` values over this `kc` range, for undoing `B`'s
+    // bias. `A_view` is already restricted to the range, so `kc` is its width.
+    int32_t a_rowsum[kNR] = {};
+    if constexpr (GEMMA_MM_I8_BIASED_B) {
+      a_rowsum[0] = A_view.RowSum(imc + 0, kc);
+      if constexpr (kRowsAC > 1) a_rowsum[1] = A_view.RowSum(imc + 1, kc);
+      if constexpr (kRowsAC > 2) a_rowsum[2] = A_view.RowSum(imc + 2, kc);
+      if constexpr (kRowsAC > 3) a_rowsum[3] = A_view.RowSum(imc + 3, kc);
+    }
+
+    MMI8StoreHorizontalSumsIntoC<kRowsAC> horz;
+    const hn::Full128<int32_t> d4i;
+    hn::Vec<decltype(d4i)> sum0, sum1, sum2, sum3;
+    horz.Reduce4x4(di32, C00, C01, C02, C03, C10, C11, C12, C13, C20, C21, C22,
+                   C23, C30, C31, C32, C33, sum0, sum1, sum2, sum3);
+    horz.Store(d4i, sum0, sum1, sum2, sum3, A_view.scale, a_rowsum, b_scale,
+               add, imc, tag, C_MC_NR);
+  }
+
+  // As `MMKernel::A2C0`.
+  template <class Tag, class CView>
+  static HWY_INLINE void A2C0(const AView A_view,
+                              const StridedView<int8_t> B_view,
+                              const float* HWY_RESTRICT b_scale, size_t mr,
+                              const IndexRange& range_mc, size_t kc,
+                              const float* HWY_RESTRICT add, Tag tag,
+                              CView C_MC_NR) {
+    HWY_DASSERT(1 <= mr && mr <= kMaxMR);
+    const size_t mc = range_mc.Num();
+    size_t imc = 0;
+
+    if (HWY_UNLIKELY(mr == 1)) {
+      for (; imc < mc; ++imc) {
+        LoopKC<1>(A_view, B_view, b_scale, imc, kc, add, tag, C_MC_NR);
+      }
+      return;
+    }
+
+    if (HWY_UNLIKELY(mr == 2)) {
+      if (HWY_LIKELY(mc >= 2)) {
+        for (; imc <= mc - 2; imc += 2) {
+          LoopKC<2>(A_view, B_view, b_scale, imc, kc, add, tag, C_MC_NR);
+        }
+      }
+      if (HWY_UNLIKELY(imc != mc)) {
+        LoopKC<1>(A_view, B_view, b_scale, imc, kc, add, tag, C_MC_NR);
+      }
+      return;
+    }
+
+    HWY_DASSERT(mr == 4);
+    if (HWY_LIKELY(mc >= 4)) {
+      for (; imc <= mc - 4; imc += 4) {
+        LoopKC<4>(A_view, B_view, b_scale, imc, kc, add, tag, C_MC_NR);
+      }
+    }
+    const size_t remainder_mc = mc - imc;
+    HWY_DASSERT(remainder_mc < 4);
+    if (HWY_UNLIKELY(remainder_mc & 2)) {
+      LoopKC<2>(A_view, B_view, b_scale, imc, kc, add, tag, C_MC_NR);
+      imc += 2;
+    }
+    if (HWY_UNLIKELY(remainder_mc & 1)) {
+      LoopKC<1>(A_view, B_view, b_scale, imc, kc, add, tag, C_MC_NR);
+      imc += 1;
+    }
+    HWY_DASSERT(imc == mc);
+  }
+};  // MMI8Kernel
+
+//------------------------------------------------------------------------------
+// Quantization
+
+// Loads one vector of F32 from F32 or BF16 `A`, so that quantization can read
+// activations in whichever format the caller already has.
+template <class DF, typename TA, class VF = hn::Vec<DF>>
+static HWY_INLINE VF LoadF32(DF df, const TA* HWY_RESTRICT p) {
+  if constexpr (IsF32<TA>()) {
+    return hn::LoadU(df, p);
+  } else {
+    static_assert(IsBF16<TA>());
+    return hn::PromoteTo(df, hn::LoadU(hn::Rebind<BF16, DF>(), p));
+  }
+}
+
+template <class DF, typename TA, class VF = hn::Vec<DF>>
+static HWY_INLINE VF LoadNF32(DF df, const TA* HWY_RESTRICT p, size_t n) {
+  if constexpr (IsF32<TA>()) {
+    return hn::LoadN(df, p, n);
+  } else {
+    static_assert(IsBF16<TA>());
+    return hn::PromoteTo(df, hn::LoadN(hn::Rebind<BF16, DF>(), p, n));
+  }
+}
+
+// Quantizes one row of `k` activations to symmetric int8, returning the
+// dequantization scale. Also writes `k + 1` prefix sums of the quantized
+// values (when `B` is biased), which the kernel uses to undo that bias for
+// whichever `kc` range it is working on. `out` is zero-padded to `padded_k`.
+template <typename TA>
+static HWY_INLINE float QuantizeRowA(const TA* HWY_RESTRICT in, size_t k,
+                                     MMI8AT* HWY_RESTRICT out,
+                                     int32_t* HWY_RESTRICT prefix,
+                                     size_t padded_k) {
+  const hn::ScalableTag<float> df;
+  const hn::Rebind<int32_t, decltype(df)> di32;
+  const hn::Rebind<MMI8AT, decltype(df)> d8;
+  using VF = hn::Vec<decltype(df)>;
+  const size_t NF = hn::Lanes(df);
+
+  VF vmax = hn::Zero(df);
+  size_t i = 0;
+  if (k >= NF) {
+    for (; i <= k - NF; i += NF) {
+      vmax = hn::Max(vmax, hn::Abs(LoadF32(df, in + i)));
+    }
+  }
+  if (i != k) {
+    vmax = hn::Max(vmax, hn::Abs(LoadNF32(df, in + i, k - i)));
+  }
+  const float amax = hn::ReduceMax(df, vmax);
+
+  const float scale = (amax == 0.0f) ? 1.0f : amax / kMMI8Max;
+  const float inv_scale = (amax == 0.0f) ? 0.0f : kMMI8Max / amax;
+  const VF vinv = hn::Set(df, inv_scale);
+
+  i = 0;
+  if (k >= NF) {
+    for (; i <= k - NF; i += NF) {
+      const auto q = hn::NearestInt(hn::Mul(LoadF32(df, in + i), vinv));
+      hn::StoreU(hn::DemoteTo(d8, q), d8, out + i);
+    }
+  }
+  for (; i < k; ++i) {
+    const float in_f = hwy::ConvertScalarTo<float>(in[i]);
+    out[i] = static_cast<MMI8AT>(std::lroundf(in_f * inv_scale));
+  }
+  for (; i < padded_k; ++i) {
+    out[i] = static_cast<MMI8AT>(0);
+  }
+
+  if constexpr (GEMMA_MM_I8_BIASED_B) {
+    // Scalar, but only `M * K` additions per MatMul, i.e. the same order as
+    // the quantization itself and negligible next to `M * K * N` products.
+    int32_t sum = 0;
+    prefix[0] = 0;
+    for (size_t j = 0; j < k; ++j) {
+      sum += out[j];
+      prefix[j + 1] = sum;
+    }
+  }
+  return scale;
+}
+
+// Storage for quantized `A`, reused across `MatMulI8` calls. Analogous to
+// `MMEntireA`, but sized by the caller because this is a prototype and
+// `MatMulEnv` does not know about int8 yet.
+class MMI8AStorage {
+ public:
+  // `prefix_` is `K + 1` per row, which is simple but the largest cost here.
+  // Production would instead compute one sum per (row, kc range) once the
+  // config is known, which is `NumTasks()` rather than `K` per row.
+  MMI8AStorage(size_t max_M, size_t max_K, const Allocator& allocator)
+      : data_("A_i8", Extents2D(max_M, max_K), allocator, MatPadding::kOdd),
+        prefix_stride_(hwy::RoundUpTo(max_K + 1, HWY_ALIGNMENT / 4)),
+        prefix_((GEMMA_MM_I8_BIASED_B ? max_M : 1) * prefix_stride_),
+        scale_(max_M) {}
+
+  MMI8AView View(const Extents2D& extents) {
+    HWY_DASSERT(extents.rows <= data_.Rows());
+    HWY_DASSERT(extents.cols <= data_.Cols());
+    return MMI8AView{
+        StridedView<MMI8AT>(HWY_RCAST_ALIGNED(MMI8AT*, data_.Row(0)),
+                            extents.cols, data_.Stride()),
+        scale_.data(), prefix_.data(), prefix_stride_};
+  }
+
+  float* HWY_RESTRICT scale() { return scale_.data(); }
+  int32_t* HWY_RESTRICT prefix(size_t row) {
+    return prefix_.data() + (GEMMA_MM_I8_BIASED_B ? row : 0) * prefix_stride_;
+  }
+  size_t Stride() const { return data_.Stride(); }
+
+ private:
+  MatStorageT<uint8_t> data_;
+  size_t prefix_stride_;
+  hwy::AlignedVector<int32_t> prefix_;
+  hwy::AlignedVector<float> scale_;
+};
+
+// Quantizes all `M x K` of `A` into `storage`, in parallel over rows.
+// This replaces `MMDecompress::DecompressA` and is the same order of cost:
+// one pass over `A`, once per `MatMul` rather than per B tile.
+template <typename TA>
+static HWY_NOINLINE MMI8AView QuantizeA(const MatPtrT<TA>& A,
+                                        MMI8AStorage& storage,
+                                        ThreadingContext& ctx,
+                                        size_t cluster_idx) {
+  const MMI8AView view = storage.View(A.Extents());
+  const size_t k = A.Cols();
+  const size_t padded_k = hwy::RoundUpTo(k, hn::Lanes(hn::ScalableTag<int8_t>()));
+  float* HWY_RESTRICT scale = storage.scale();
+  const float a_scale = A.Scale();
+
+  ParallelFor(Parallelism::kFlat, A.Rows(), ctx, cluster_idx,
+              Callers::kMMQuantizeA,
+              [&](size_t r, size_t /*worker*/) HWY_ATTR {
+                if (MMI8CanRotate(k)) {
+                  hwy::AlignedVector<float> rotated(padded_k);
+                  for (size_t c = 0; c < k; ++c) {
+                    rotated[c] = hwy::ConvertScalarTo<float>(A.Row(r)[c]);
+                  }
+                  MMI8Rotate(rotated.data(), k);
+                  scale[r] =
+                      a_scale * QuantizeRowA(rotated.data(), k,
+                                             view.data.Row(r), storage.prefix(r),
+                                             padded_k);
+                } else {
+                  scale[r] =
+                      a_scale * QuantizeRowA(A.Row(r), k, view.data.Row(r),
+                                             storage.prefix(r), padded_k);
+                }
+              });
+  return view;
+}
+
+// Symmetric int8 quantization of already-transposed `B`, i.e. `N` rows of `K`.
+// Fills `data` (biased by 128 if `GEMMA_MM_I8_BIASED_B`, zero-padded to its
+// stride) and `scale`. Called once per weight matrix, so not performance-
+// critical.
+static HWY_NOINLINE MMI8B PackB(const MatPtrT<float>& B_f32,
+                                MatPtrT<int8_t>& data,
+                                float* HWY_RESTRICT scale,
+                                ThreadingContext& ctx) {
+  const size_t k = B_f32.Cols();
+  const float b_scale = B_f32.Scale();
+
+  ParallelFor(Parallelism::kFlat, B_f32.Rows(), ctx, /*cluster_idx=*/0,
+              Callers::kTest, [&](size_t r, size_t /*worker*/) HWY_ATTR {
+                const float* HWY_RESTRICT in = B_f32.Row(r);
+                hwy::AlignedVector<float> rotated;
+                if (MMI8CanRotate(k)) {
+                  rotated.resize(k);
+                  hwy::CopyBytes(in, rotated.data(), k * sizeof(float));
+                  MMI8Rotate(rotated.data(), k);
+                  in = rotated.data();
+                }
+                float amax = 0.0f;
+                for (size_t c = 0; c < k; ++c) {
+                  amax = HWY_MAX(amax, hwy::ScalarAbs(in[c]));
+                }
+                const float s = (amax == 0.0f) ? 1.0f : amax / kMMI8Max;
+                const float inv = (amax == 0.0f) ? 0.0f : kMMI8Max / amax;
+                MMI8BT* HWY_RESTRICT out =
+                    HWY_RCAST_ALIGNED(MMI8BT*, data.Row(r));
+                for (size_t c = 0; c < k; ++c) {
+                  const int32_t q =
+                      static_cast<int32_t>(std::lroundf(in[c] * inv));
+                  HWY_DASSERT(-127 <= q && q <= 127);
+                  out[c] = static_cast<MMI8BT>(
+                      q + (GEMMA_MM_I8_BIASED_B ? 128 : 0));
+                }
+                for (size_t c = k; c < data.Stride(); ++c) {
+                  out[c] = static_cast<MMI8BT>(0);
+                }
+                scale[r] = b_scale * s;
+              });
+
+  return MMI8B{&data, scale};
+}
+
+//------------------------------------------------------------------------------
+// Entry point
+
+// As `MatMul`, but `A` is quantized on the fly and `B` was packed by `PackB`.
+// Reuses the same blocking, parallelization and autotuning as `MatMul`; only
+// the kernel and operand types differ. `env` must not be shared with
+// (BF16) `MatMul` calls of the same shape, because the autotuner is keyed on
+// shape alone and the two kernels prefer different configs.
+template <typename TA, typename TC>
+HWY_NOINLINE MMPerKey* MatMulI8(const MatPtrT<TA>& A, const MMI8B& B,
+                                const float* HWY_RESTRICT add, MatMulEnv& env,
+                                MatPtrT<TC>& C, MMI8AStorage& a_storage,
+                                MMOptions options = MMOptions()) {
+  const size_t cluster_idx = options.cluster_idx;
+  HWY_DASSERT(cluster_idx < env.row_ptrs.size());
+  GCPP_ZONE(env.ctx, env.ctx.Worker(cluster_idx), Zones::kMMMatMul);
+
+  RowPtrs<TC> C_rows = GetOrSetTempRowPtrs(C, env.row_ptrs[cluster_idx]);
+
+  const size_t M = A.Rows();
+  const size_t K = A.Cols();
+  const size_t N = B.Rows();
+  const size_t num_B = 1;
+
+  const CacheInfo& cache = env.ctx.cache_info;
+  MMPerKey& per_key = MMImpl::FindOrAddPerKey(
+      M, K, N, num_B, cache.VectorBytes(), env.per_cluster[cluster_idx]);
+
+  // Outside the timed section, as `MMDecompress::MaybeDecompressA`.
+  const MMI8AView A_view = QuantizeA(A, a_storage, env.ctx, cluster_idx);
+
+  const MMI8B* B2 = nullptr;  // required for type matching
+
+  // Scales are per row/column, hence folded into `A_view.scale` and
+  // `B.scale`; the scalar `MMArgs::scale_A` is unused.
+  MMAutoTune<MMConfig>& tuner = per_key.autotune;
+  if (HWY_LIKELY(tuner.Best())) {
+    const MMArgs args(env, M, K, N, /*scale_A=*/1.0f, add, options, tuner,
+                      *tuner.Best());
+    MMLoops::Dispatch<MMI8Kernel>(A_view, B, B2, C_rows, args);
+    return &per_key;
+  }
+
+  if (HWY_UNLIKELY(!tuner.HasCandidates())) {
+    HWY_ASSERT(K == B.Cols());
+    HWY_ASSERT(M <= kMaxBatchSize);
+    HWY_ASSERT(N % kNR == 0);
+    tuner.SetCandidates(
+        MMCandidates(cache, M, K, N, num_B, sizeof(TC), env.print_config));
+  }
+
+  const MMConfig& cfg = tuner.NextConfig();
+  const MMArgs args(env, M, K, N, /*scale_A=*/1.0f, add, options, tuner, cfg);
+
+  const uint64_t t0 = hwy::timer::Start();
+  MMLoops::Dispatch<MMI8Kernel>(A_view, B, B2, C_rows, args);
+  MMImpl::NotifyAutotuneResult(env, M, K, N, num_B, t0, tuner, cfg);
+
+  return &per_key;
+}
+
+// As `TwoMatMul`: computes `A * B1` into `C` and `A * B2` into a per-worker
+// tile, passing both to `options.func`. Used by gated FFNs.
+static HWY_NOINLINE MMPerKey* TwoMatMulI8(const MatPtrT<BF16>& A,
+                                          const MMI8B& B1,
+                                          const MMI8B& B2, MatMulEnv& env,
+                                          MatPtrT<BF16>& C,
+                                          MMI8AStorage& a_storage,
+                                          MMOptions options) {
+  const size_t cluster_idx = options.cluster_idx;
+  HWY_DASSERT(cluster_idx < env.row_ptrs.size());
+  GCPP_ZONE(env.ctx, env.ctx.Worker(cluster_idx), Zones::kMMTwoMatMul);
+  HWY_DASSERT(options.func != nullptr);  // no other way to get access to C2.
+
+  RowPtrs<BF16> C_rows = GetOrSetTempRowPtrs(C, env.row_ptrs[cluster_idx]);
+
+  const size_t M = A.Rows();
+  const size_t K = A.Cols();
+  const size_t N = B1.Rows();
+  const size_t num_B = 2;
+
+  const CacheInfo& cache = env.ctx.cache_info;
+  MMPerKey& per_key = MMImpl::FindOrAddPerKey(
+      M, K, N, num_B, cache.VectorBytes(), env.per_cluster[cluster_idx]);
+
+  const MMI8AView A_view = QuantizeA(A, a_storage, env.ctx, cluster_idx);
+
+  MMAutoTune<MMConfig>& tuner = per_key.autotune;
+  if (HWY_LIKELY(tuner.Best())) {
+    const MMArgs args(env, M, K, N, /*scale_A=*/1.0f, /*add=*/nullptr, options,
+                      tuner, *tuner.Best());
+    MMLoops::Dispatch<MMI8Kernel>(A_view, B1, &B2, C_rows, args);
+    return &per_key;
+  }
+
+  if (HWY_UNLIKELY(!tuner.HasCandidates())) {
+    HWY_ASSERT(K == B1.Cols());
+    HWY_ASSERT(K == B2.Cols());
+    HWY_ASSERT(M <= kMaxBatchSize);
+    HWY_ASSERT(N % kNR == 0);
+    const size_t max_M = MMKeys::BucketM(M);
+    tuner.SetCandidates(MMCandidates(cache, max_M, K, N, num_B, sizeof(BF16),
+                                     env.print_config));
+  }
+
+  const MMConfig& cfg = tuner.NextConfig();
+  const MMArgs args(env, M, K, N, /*scale_A=*/1.0f, /*add=*/nullptr, options,
+                    tuner, cfg);
+  const uint64_t t0 = hwy::timer::Start();
+  MMLoops::Dispatch<MMI8Kernel>(A_view, B1, &B2, C_rows, args);
+  MMImpl::NotifyAutotuneResult(env, M, K, N, num_B, t0, tuner, cfg);
+
+  return &per_key;
+}
+
+// NOLINTNEXTLINE(google-readability-namespace-comments)
+}  // namespace HWY_NAMESPACE
+}  // namespace gcpp
+HWY_AFTER_NAMESPACE();
+
+#endif  // NOLINT
diff --git a/ops/matmul_i8_model-inl.h b/ops/matmul_i8_model-inl.h
new file mode 100644
index 00000000..736ac4c6
--- /dev/null
+++ b/ops/matmul_i8_model-inl.h
@@ -0,0 +1,302 @@
+// Copyright 2025 Google LLC
+// SPDX-License-Identifier: Apache-2.0
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Experiment harness that routes the model's MatMuls through the W8A8 kernel
+// in `ops/matmul_i8-inl.h`, to measure end-to-end quality. Weights are
+// quantized lazily on first use and cached, keyed by their data pointer, so
+// this needs no changes to the loading path.
+//
+// NOT a production integration:
+//  - quantizing from whatever the file holds (e.g. SFP) stacks a second
+//    quantization on top of the first; a real path would quantize the
+//    original checkpoint;
+//  - the int8 and BF16 kernels share `MatMulEnv`'s autotune keys, which are
+//    shape-only, so a model that mixes them will mis-tune;
+//  - the cache is a process-wide singleton and never freed.
+//
+// Enabled by environment variables, so no CLI plumbing is needed:
+//   GEMMA_MM_I8=1            route eligible MatMuls through the int8 kernel
+//   GEMMA_MM_I8_MIN_K=<n>    leave tensors with K < n in their original format
+//   GEMMA_MM_I8_SKIP_ROWS=<n> leave tensors with N >= n alone (e.g. the vocab-
+//                            sized logits projection, the usual first thing to
+//                            exclude from W8A8)
+//   GEMMA_MM_I8_INCLUDE=<list> only quantize tensor names containing one of the
+//                            comma-separated substrings
+//   GEMMA_MM_I8_EXCLUDE=<list> leave matching tensor names in their old format
+//   GEMMA_MM_I8_ROTATE=1       apply matching block-Hadamard rotations to A/B
+//   GEMMA_MM_I8_VERBOSE=1    log each tensor as it is quantized
+
+#include <stddef.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <memory>
+#include <mutex>  // NOLINT
+#include <unordered_map>
+
+#include "ops/matmul.h"
+#include "util/mat.h"
+#include "util/threading_context.h"
+#include "hwy/base.h"
+
+// Include guard for (potentially) SIMD code.
+#if defined(THIRD_PARTY_GEMMA_CPP_MATMUL_I8_MODEL_TOGGLE) == \
+    defined(HWY_TARGET_TOGGLE)
+#ifdef THIRD_PARTY_GEMMA_CPP_MATMUL_I8_MODEL_TOGGLE
+#undef THIRD_PARTY_GEMMA_CPP_MATMUL_I8_MODEL_TOGGLE
+#else
+#define THIRD_PARTY_GEMMA_CPP_MATMUL_I8_MODEL_TOGGLE
+#endif
+
+#include "hwy/highway.h"
+// After highway.h
+#include "compression/compress-inl.h"
+#include "ops/matmul_i8-inl.h"
+
+HWY_BEFORE_NAMESPACE();
+namespace gcpp {
+namespace HWY_NAMESPACE {
+namespace hn = hwy::HWY_NAMESPACE;
+
+// Reads an integer environment variable, or returns `fallback`.
+static inline size_t MMI8EnvSize(const char* name, size_t fallback) {
+  const char* s = getenv(name);
+  if (s == nullptr || *s == '\0') return fallback;
+  const long long v = atoll(s);  // NOLINT
+  return v < 0 ? fallback : static_cast<size_t>(v);
+}
+
+// True if `name` contains any non-empty comma-separated token in `list`.
+static inline bool MMI8NameMatches(const char* name, const char* list) {
+  if (list == nullptr || *list == '\0') return false;
+  for (const char* begin = list; *begin != '\0';) {
+    const char* end = strchr(begin, ',');
+    if (end == nullptr) end = begin + strlen(begin);
+    const size_t len = static_cast<size_t>(end - begin);
+    if (len != 0) {
+      for (const char* at = name; *at != '\0'; ++at) {
+        if (strncmp(at, begin, len) == 0) return true;
+      }
+    }
+    begin = *end == '\0' ? end : end + 1;
+  }
+  return false;
+}
+
+// Quantizes one row of `k` floats to symmetric int8, biased by 128 if
+// `GEMMA_MM_I8_BIASED_B`. Returns the dequantization scale.
+static HWY_INLINE float PackBRow(const float* HWY_RESTRICT in, size_t k,
+                                 MMI8BT* HWY_RESTRICT out, size_t padded_k) {
+  const hn::ScalableTag<float> df;
+  const hn::Rebind<int32_t, decltype(df)> di32;
+  const hn::Rebind<int8_t, decltype(df)> di8;
+  using VF = hn::Vec<decltype(df)>;
+  const size_t NF = hn::Lanes(df);
+
+  VF vmax = hn::Zero(df);
+  size_t i = 0;
+  if (k >= NF) {
+    for (; i <= k - NF; i += NF) {
+      vmax = hn::Max(vmax, hn::Abs(hn::LoadU(df, in + i)));
+    }
+  }
+  if (i != k) vmax = hn::Max(vmax, hn::Abs(hn::LoadN(df, in + i, k - i)));
+  const float amax = hn::ReduceMax(df, vmax);
+
+  const float scale = (amax == 0.0f) ? 1.0f : amax / kMMI8Max;
+  const float inv = (amax == 0.0f) ? 0.0f : kMMI8Max / amax;
+  const VF vinv = hn::Set(df, inv);
+  // Store as int8 and add the bias afterwards: `DemoteTo` to u8 would saturate
+  // negative values to zero.
+  const auto vbias = hn::Set(di32, GEMMA_MM_I8_BIASED_B ? 128 : 0);
+
+  i = 0;
+  if (k >= NF) {
+    for (; i <= k - NF; i += NF) {
+      const auto q = hn::NearestInt(hn::Mul(hn::LoadU(df, in + i), vinv));
+      // Bias in the int32 domain, then narrow; `DemoteTo` saturates, and
+      // `q + 128` is within [1, 255] so nothing is clamped.
+      if constexpr (GEMMA_MM_I8_BIASED_B) {
+        const hn::Rebind<uint8_t, decltype(df)> du8;
+        hn::StoreU(hn::DemoteTo(du8, hn::Add(q, vbias)), du8,
+                   HWY_RCAST_ALIGNED(uint8_t*, out) + i);
+      } else {
+        hn::StoreU(hn::DemoteTo(di8, q), di8,
+                   HWY_RCAST_ALIGNED(int8_t*, out) + i);
+      }
+    }
+  }
+  for (; i < k; ++i) {
+    const int32_t q = static_cast<int32_t>(std::lroundf(in[i] * inv));
+    out[i] = static_cast<MMI8BT>(q + (GEMMA_MM_I8_BIASED_B ? 128 : 0));
+  }
+  for (; i < padded_k; ++i) out[i] = static_cast<MMI8BT>(0);
+  return scale;
+}
+
+// Process-wide cache of int8 weights, keyed by the tensor's data pointer.
+class MMI8WeightCache {
+ public:
+  static MMI8WeightCache& Get() {
+    static MMI8WeightCache cache;
+    return cache;
+  }
+
+  bool Enabled() const { return enabled_; }
+
+  // Returns the int8 form of `B`, quantizing and caching on first use, or
+  // nullptr if this tensor is not eligible (see the environment variables).
+  template <typename TB>
+  const MMI8B* Lookup(const MatPtrT<TB>& B, ThreadingContext& ctx) {
+    const size_t N = B.Rows();
+    const size_t K = B.Cols();
+    if (K < min_k_ || N >= skip_rows_ || (N % kNR) != 0) return nullptr;
+    if (include_ != nullptr && *include_ != '\0' &&
+        !MMI8NameMatches(B.Name(), include_)) {
+      return nullptr;
+    }
+    if (MMI8NameMatches(B.Name(), exclude_)) return nullptr;
+
+    const void* key = B.RowBytes(0);
+    std::lock_guard<std::mutex> lock(mutex_);
+    auto it = map_.find(key);
+    if (it != map_.end()) return it->second ? &it->second->b : nullptr;
+
+    auto entry = std::make_unique<Entry>(B, ctx.allocator);
+    Quantize(B, *entry);
+    if (verbose_) {
+      fprintf(stderr, "MM.I8: quantized %-16s %6zu x %6zu\n", B.Name(), N, K);
+    }
+    const MMI8B* result = &entry->b;
+    map_[key] = std::move(entry);
+    return result;
+  }
+
+  // Storage for the quantized `A`, grown on demand. `MatMul` for a given
+  // `MatMulEnv` is not called concurrently, and this experiment runs a single
+  // cluster, so one instance suffices.
+  MMI8AStorage& AStorage(size_t M, size_t K, const Allocator& allocator) {
+    if (a_ == nullptr || M > a_max_M_ || K > a_max_K_) {
+      a_max_M_ = HWY_MAX(M, a_max_M_);
+      a_max_K_ = HWY_MAX(K, a_max_K_);
+      a_ = std::make_unique<MMI8AStorage>(a_max_M_, a_max_K_, allocator);
+    }
+    return *a_;
+  }
+
+ private:
+  struct Entry {
+    Entry(const MatPtr& B, const Allocator& allocator)
+        : data("B_i8", Extents2D(B.Rows(), B.Cols()), allocator,
+               MatPadding::kOdd),
+          scale(B.Rows()) {
+      b = MMI8B{&data, scale.data()};
+    }
+    MatStorageT<int8_t> data;
+    hwy::AlignedVector<float> scale;
+    MMI8B b;
+  };
+
+  MMI8WeightCache()
+      : enabled_(MMI8EnvSize("GEMMA_MM_I8", 0) != 0),
+        verbose_(MMI8EnvSize("GEMMA_MM_I8_VERBOSE", 0) != 0),
+        min_k_(MMI8EnvSize("GEMMA_MM_I8_MIN_K", 0)),
+        skip_rows_(MMI8EnvSize("GEMMA_MM_I8_SKIP_ROWS", ~size_t{0})),
+        include_(getenv("GEMMA_MM_I8_INCLUDE")),
+        exclude_(getenv("GEMMA_MM_I8_EXCLUDE")) {}
+
+  // Serial (the caller may already be inside a parallel region), but
+  // vectorized, so a 2B-parameter model takes a few seconds in total.
+  template <typename TB>
+  void Quantize(const MatPtrT<TB>& B, Entry& entry) {
+    const hn::ScalableTag<float> df;
+    const size_t K = B.Cols();
+    const size_t padded_k = hwy::RoundUpTo(K, hn::Lanes(df));
+    hwy::AlignedVector<float> row(padded_k + hn::Lanes(df));
+    const PackedSpan<const TB> span = B.PaddedSpan();
+    const float b_scale = B.Scale();
+
+    for (size_t r = 0; r < B.Rows(); ++r) {
+      DecompressAndZeroPad(df, span, r * B.Stride(), row.data(), K);
+      if (MMI8CanRotate(K)) MMI8Rotate(row.data(), K);
+      MMI8BT* HWY_RESTRICT out =
+          HWY_RCAST_ALIGNED(MMI8BT*, entry.data.Row(r));
+      entry.scale[r] =
+          b_scale * PackBRow(row.data(), K, out, entry.data.Stride());
+    }
+  }
+
+  bool enabled_;
+  bool verbose_;
+  size_t min_k_;
+  size_t skip_rows_;
+  const char* include_;
+  const char* exclude_;
+
+  std::mutex mutex_;
+  std::unordered_map<const void*, std::unique_ptr<Entry>> map_;
+
+  std::unique_ptr<MMI8AStorage> a_;
+  size_t a_max_M_ = 0;
+  size_t a_max_K_ = 0;
+};
+
+// As `MaybeMatMulI8`, for the fused gated-FFN pair. Both operands must be
+// eligible, else we fall back so that the pair stays consistent.
+static inline MMPerKey* MaybeTwoMatMulI8(const MatPtrT<BF16>& A,
+                                         const MatPtr& B1, const MatPtr& B2,
+                                         MatMulEnv& env, MatPtrT<BF16>& C,
+                                         const MMOptions& options) {
+  MMI8WeightCache& cache = MMI8WeightCache::Get();
+  if (!cache.Enabled()) return nullptr;
+  return CallUpcastedSame(
+      &B1, &B2, [&](const auto* B1_t, const auto* B2_t) -> MMPerKey* {
+        const MMI8B* i8_1 = cache.Lookup(*B1_t, env.ctx);
+        if (i8_1 == nullptr) return nullptr;
+        const MMI8B* i8_2 = cache.Lookup(*B2_t, env.ctx);
+        if (i8_2 == nullptr) return nullptr;
+        MMI8AStorage& a_storage =
+            cache.AStorage(A.Rows(), A.Cols(), env.ctx.allocator);
+        return TwoMatMulI8(A, *i8_1, *i8_2, env, C, a_storage, options);
+      });
+}
+
+// If the int8 path is enabled and `B` is eligible, computes `C = A * B + add`
+// with the W8A8 kernel and returns its autotune state; else returns nullptr so
+// the caller falls back to `MatMulStatic`.
+template <typename TA, typename TB, typename TC>
+MMPerKey* MaybeMatMulI8(const MatPtrT<TA>& A, const MatPtrT<TB>& B,
+                        const float* HWY_RESTRICT add, MatMulEnv& env,
+                        MatPtrT<TC>& C, const MMOptions& options) {
+  MMI8WeightCache& cache = MMI8WeightCache::Get();
+  if (!cache.Enabled()) return nullptr;
+  // `TwoMatMul`'s fused second output is not wired up here.
+  if (options.func != nullptr) return nullptr;
+  const MMI8B* B_i8 = cache.Lookup(B, env.ctx);
+  if (B_i8 == nullptr) return nullptr;
+
+  MMI8AStorage& a_storage =
+      cache.AStorage(A.Rows(), A.Cols(), env.ctx.allocator);
+  return MatMulI8(A, *B_i8, add, env, C, a_storage, options);
+}
+
+// NOLINTNEXTLINE(google-readability-namespace-comments)
+}  // namespace HWY_NAMESPACE
+}  // namespace gcpp
+HWY_AFTER_NAMESPACE();
+
+#endif  // NOLINT
diff --git a/ops/matmul_i8_test.cc b/ops/matmul_i8_test.cc
new file mode 100644
index 00000000..1c8be2ff
--- /dev/null
+++ b/ops/matmul_i8_test.cc
@@ -0,0 +1,334 @@
+// Copyright 2025 Google LLC
+// SPDX-License-Identifier: Apache-2.0
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//      http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Correctness of the W8A8 kernel in `ops/matmul_i8-inl.h`. The reference is
+// computed in F64 from the *quantized* operands, so this checks the kernel's
+// arithmetic (accumulation, remainder handling, the u8 bias correction, and
+// the MMSetC/MMAddC split across kc ranges) rather than quantization error.
+//
+// Built twice, with `GEMMA_MM_I8_FORCE_BIASED_B` 0 and 1, so that the x86
+// biased-u8 path is covered on non-x86 hosts too.
+
+#include <stddef.h>
+#include <stdint.h>
+#include <stdio.h>
+
+#include <cmath>
+#include <vector>
+
+#include "compression/types.h"  // GEMMA_DISABLED_TARGETS
+#ifndef HWY_DISABLED_TARGETS
+#define HWY_DISABLED_TARGETS GEMMA_DISABLED_TARGETS
+#endif  // HWY_DISABLED_TARGETS
+
+#include "ops/matmul.h"
+#include "util/basics.h"
+#include "util/mat.h"
+#include "util/threading_context.h"
+#include "hwy/aligned_allocator.h"
+
+// clang-format off
+#undef HWY_TARGET_INCLUDE
+#define HWY_TARGET_INCLUDE "ops/matmul_i8_test.cc"  // NOLINT
+// clang-format on
+#include "hwy/foreach_target.h"  // IWYU pragma: keep
+#include "hwy/highway.h"
+// After highway.h
+#include "compression/compress-inl.h"
+#include "ops/matmul-inl.h"
+#include "ops/matmul_i8-inl.h"
+
+HWY_BEFORE_NAMESPACE();
+namespace gcpp {
+
+// Not in HWY_NAMESPACE: the `HWY_ONCE` section below must read the same
+// instance that the dispatched target wrote to.
+extern size_t g_failures;
+
+namespace HWY_NAMESPACE {
+namespace hn = hwy::HWY_NAMESPACE;
+
+class Rng {
+ public:
+  explicit Rng(uint64_t seed) : state_(seed * 6364136223846793005ull + 1) {}
+  float Normal() {
+    float sum = 0.0f;
+    for (int i = 0; i < 4; ++i) sum += Uniform();
+    return (sum - 2.0f) * 1.732f;
+  }
+
+ private:
+  float Uniform() {
+    state_ = state_ * 6364136223846793005ull + 1442695040888963407ull;
+    return static_cast<float>((state_ >> 40) & 0xFFFFFF) / 16777216.0f;
+  }
+  uint64_t state_;
+};
+
+void TestRotationPreservesDotProducts() {
+  constexpr size_t kSize = 2 * kMMI8RotateBlock;
+  std::vector<float> a(kSize);
+  std::vector<float> b(kSize);
+  Rng rng(123);
+  double expected = 0.0;
+  for (size_t i = 0; i < kSize; ++i) {
+    a[i] = rng.Normal();
+    b[i] = rng.Normal();
+    expected += static_cast<double>(a[i]) * b[i];
+  }
+
+  MMI8Rotate(a.data(), a.size());
+  MMI8Rotate(b.data(), b.size());
+  double actual = 0.0;
+  for (size_t i = 0; i < kSize; ++i) {
+    actual += static_cast<double>(a[i]) * b[i];
+  }
+
+  const double relative =
+      hwy::ScalarAbs(actual - expected) /
+      HWY_MAX(1.0, hwy::ScalarAbs(expected));
+  if (relative > 1E-6) {
+    ++g_failures;
+    printf("FAIL rotation dot-product relative error %.3e\n", relative);
+  } else {
+    printf("  ok rotation preserves dot products (relative error %.3e)\n",
+           relative);
+  }
+}
+
+// Fills A and B. Row magnitudes deliberately vary by up to 7x, so that a
+// mixed-up per-row scale index would show up.
+void FillOperands(size_t M, size_t K, size_t N, MatStorageT<float>& A_f32,
+                  MatStorageT<float>& B_f32, float b_mean = 0.0f) {
+  Rng rng(M * 131 + K * 17 + N);
+  for (size_t r = 0; r < M; ++r) {
+    float* row = A_f32.Row(r);
+    const float row_scale = 0.01f * static_cast<float>(1 + (r % 7));
+    for (size_t c = 0; c < K; ++c) row[c] = rng.Normal() * row_scale;
+    for (size_t c = K; c < A_f32.Stride(); ++c) row[c] = 0.0f;
+  }
+  for (size_t r = 0; r < N; ++r) {
+    float* row = B_f32.Row(r);
+    const float row_scale = 0.5f * static_cast<float>(1 + (r % 5));
+    // A nonzero mean makes the per-channel sums of the quantized weights
+    // large. Correcting `B`'s bias once over the whole `K` (rather than per
+    // `kc` range) would then write intermediates to `C` that are far larger
+    // than the result, which is unrecoverable when `C` is BF16.
+    for (size_t c = 0; c < K; ++c) {
+      row[c] = (rng.Normal() + b_mean) * row_scale;
+    }
+    for (size_t c = K; c < B_f32.Stride(); ++c) row[c] = 0.0f;
+  }
+}
+
+// One `M x K x N` case. `TC` is the output type; `add` exercises the bias.
+template <typename TC>
+void TestCase(size_t M, size_t K, size_t N, bool add, ThreadingContext& ctx,
+              MatMulEnv& env, MMI8AStorage& a_i8, float b_mean = 0.0f) {
+  const Allocator& allocator = ctx.allocator;
+  MatStorageT<float> A_f32("A", Extents2D(M, K), allocator, MatPadding::kOdd);
+  MatStorageT<float> B_f32("B", Extents2D(N, K), allocator, MatPadding::kOdd);
+  FillOperands(M, K, N, A_f32, B_f32, b_mean);
+  // Non-unit tensor scales, which must be folded in by QuantizeA/PackB.
+  A_f32.SetScale(0.75f);
+  B_f32.SetScale(1.25f);
+
+  MatStorageT<int8_t> B_i8("B_i8", Extents2D(N, K), allocator,
+                           MatPadding::kOdd);
+  hwy::AlignedVector<float> b_scale(N), add_row(N);
+  const MMI8B B_packed = PackB(B_f32, B_i8, b_scale.data(), ctx);
+  for (size_t n = 0; n < N; ++n) add_row[n] = 0.125f * static_cast<float>(n % 9);
+
+  MatStorageT<TC> C("C", Extents2D(M, N), allocator, MatPadding::kOdd);
+  C.AllocateAndAttachRowPtrs(env.row_ptrs);
+  // Run until autotuning settles, then check the result produced by the best
+  // config. Otherwise every call would use a different blocking, and with a
+  // BF16 `C` the number of kc ranges changes how much precision is lost.
+  MMPerKey* per_key = nullptr;
+  for (size_t iter = 0; iter < 4096; ++iter) {
+    per_key = MatMulI8(A_f32, B_packed, add ? add_row.data() : nullptr, env, C,
+                       a_i8);
+    if (per_key->autotune.Best()) break;
+  }
+  HWY_ASSERT(per_key->autotune.Best());
+  const size_t kc = per_key->autotune.Best()->KC();
+  const size_t k_ranges = per_key->autotune.Best()->RangesOfKC(K).NumTasks();
+  MatMulI8(A_f32, B_packed, add ? add_row.data() : nullptr, env, C, a_i8);
+
+  // Reference from the quantized operands. `QuantizeA` has already written
+  // them, so read them back rather than re-deriving.
+  const MMI8AView A_q = a_i8.View(Extents2D(M, K));
+  double max_abs_err = 0.0;
+  double sum_sq = 0.0;
+  for (size_t m = 0; m < M; ++m) {
+    const MMI8AT* qa = A_q.data.Row(m);
+    for (size_t n = 0; n < N; ++n) {
+      const MMI8BT* qb = HWY_RCAST_ALIGNED(const MMI8BT*, B_i8.Row(n));
+      int64_t dot = 0;
+      for (size_t k = 0; k < K; ++k) {
+        const int32_t b =
+            static_cast<int32_t>(qb[k]) - (GEMMA_MM_I8_BIASED_B ? 128 : 0);
+        dot += static_cast<int64_t>(qa[k]) * static_cast<int64_t>(b);
+      }
+      const double expected = static_cast<double>(A_q.scale[m]) * b_scale[n] *
+                                  static_cast<double>(dot) +
+                              (add ? add_row[n] : 0.0f);
+      const double actual = hwy::ConvertScalarTo<double>(C.Row(m)[n]);
+      max_abs_err = HWY_MAX(max_abs_err, hwy::ScalarAbs(actual - expected));
+      sum_sq += expected * expected;
+    }
+  }
+  // Individual outputs are sums of `K` signed products and can cancel to near
+  // zero, where an elementwise relative error is meaningless. Normalize the
+  // worst absolute error by the RMS of the expected outputs instead.
+  const double rms = std::sqrt(sum_sq / static_cast<double>(M * N));
+  const double err = (rms == 0.0) ? 0.0 : max_abs_err / rms;
+
+  // BF16 output has 8 mantissa bits, and `MMAddC` rounds once per kc range.
+  const double tolerance = IsBF16<TC>() ? 6E-2 : 1E-5;
+  const bool ok = err <= tolerance;
+  if (!ok) ++g_failures;
+  printf(
+      "%s M=%4zu K=%5zu N=%5zu add=%d TC=%-5s biasedB=%d kc=%5zu(x%zu) "
+      "err/rms=%.2e\n",
+      ok ? "  ok  " : "FAILED", M, K, N, add, TypeName<TC>(),
+      GEMMA_MM_I8_BIASED_B, kc, k_ranges, err);
+}
+
+// Control: how much precision the *existing* BF16 kernel loses when `TC` is
+// BF16 and `K` spans several kc ranges, so that `MMAddC` accumulates through
+// BF16. Reported as a reference point for the int8 kernel's BF16-output
+// tolerance, since both inherit this from `MMStoreHorizontalSumsIntoC`.
+void ControlBF16OutputError(size_t M, size_t K, size_t N,
+                            ThreadingContext& ctx, MatMulEnv& env) {
+  const Allocator& allocator = ctx.allocator;
+  MatStorageT<float> A_f32("A", Extents2D(M, K), allocator, MatPadding::kOdd);
+  MatStorageT<float> B_f32("B", Extents2D(N, K), allocator, MatPadding::kOdd);
+  FillOperands(M, K, N, A_f32, B_f32);
+
+  MatStorageT<BF16> A_bf("A_bf", Extents2D(M, K), allocator, MatPadding::kOdd);
+  MatStorageT<BF16> B_bf("B_bf", Extents2D(N, K), allocator, MatPadding::kOdd);
+  CompressWorkingSet ws;
+  ws.tls.resize(ctx.pools.MaxWorkers());
+  for (size_t r = 0; r < M; ++r) {
+    Compress(A_f32.Row(r), K, ws.tls[0], MakeSpan(A_bf.Row(r), K), 0);
+  }
+  for (size_t r = 0; r < N; ++r) {
+    Compress(B_f32.Row(r), K, ws.tls[0], MakeSpan(B_bf.Row(r), K), 0);
+  }
+
+  MatStorageT<float> C_f32("Cf", Extents2D(M, N), allocator, MatPadding::kOdd);
+  MatStorageT<BF16> C_bf("Cb", Extents2D(M, N), allocator, MatPadding::kOdd);
+  C_f32.AllocateAndAttachRowPtrs(env.row_ptrs);
+  for (size_t iter = 0; iter < 4096; ++iter) {
+    if (MatMul(A_bf, B_bf, nullptr, env, C_f32)->autotune.Best()) break;
+  }
+  MatMul(A_bf, B_bf, nullptr, env, C_f32);
+  C_bf.AllocateAndAttachRowPtrs(env.row_ptrs);
+  MMPerKey* per_key = nullptr;
+  for (size_t iter = 0; iter < 4096; ++iter) {
+    per_key = MatMul(A_bf, B_bf, nullptr, env, C_bf);
+    if (per_key->autotune.Best()) break;
+  }
+  HWY_ASSERT(per_key->autotune.Best());
+  const size_t kc = per_key->autotune.Best()->KC();
+  const size_t k_ranges = per_key->autotune.Best()->RangesOfKC(K).NumTasks();
+  MatMul(A_bf, B_bf, nullptr, env, C_bf);
+
+  double max_abs = 0.0, sum_sq = 0.0;
+  for (size_t m = 0; m < M; ++m) {
+    for (size_t n = 0; n < N; ++n) {
+      const double f = C_f32.Row(m)[n];
+      const double b = hwy::ConvertScalarTo<double>(C_bf.Row(m)[n]);
+      max_abs = HWY_MAX(max_abs, hwy::ScalarAbs(f - b));
+      sum_sq += f * f;
+    }
+  }
+  const double rms = std::sqrt(sum_sq / static_cast<double>(M * N));
+  printf(
+      "control  M=%4zu K=%5zu N=%5zu kc=%5zu(x%zu) bf16 kernel, TC=bf16 vs "
+      "TC=f32: err/rms=%.2e\n",
+      M, K, N, kc, k_ranges, rms == 0.0 ? 0.0 : max_abs / rms);
+}
+
+void TestAll() {
+  ThreadingArgs threading_args;
+  ThreadingContext ctx(threading_args);
+  MatMulEnv env(ctx);
+  printf("target=%s biasedB=%d vector bytes=%zu\n", hwy::TargetName(HWY_TARGET),
+         GEMMA_MM_I8_BIASED_B, hn::Lanes(hn::ScalableTag<uint8_t>()));
+  TestRotationPreservesDotProducts();
+
+  // `kMaxKC` is 6 KiB, so K = 20000 forces several kc ranges and thus the
+  // MMSetC-then-MMAddC path where the bias correction must be applied once.
+  MMI8AStorage a_i8(/*max_M=*/64, /*max_K=*/20096, ctx.allocator);
+
+  // Vector-length remainders: K deliberately not a multiple of 16/32/64.
+  for (size_t K : {size_t{4}, size_t{15}, size_t{16}, size_t{17}, size_t{63},
+                   size_t{64}, size_t{65}, size_t{127}, size_t{257}}) {
+    TestCase<float>(4, K, 8, /*add=*/false, ctx, env, a_i8);
+  }
+
+  // `kRowsAC` 1/2/4 and the M remainder handling in `A2C0`.
+  for (size_t M : {size_t{1}, size_t{2}, size_t{3}, size_t{4}, size_t{5},
+                   size_t{7}, size_t{8}, size_t{13}, size_t{64}}) {
+    TestCase<float>(M, 1153, 12, /*add=*/true, ctx, env, a_i8);
+  }
+
+  // N is required to be a multiple of kNR.
+  for (size_t N : {size_t{4}, size_t{8}, size_t{16}, size_t{100},
+                   size_t{1536}}) {
+    TestCase<float>(4, 512, N, /*add=*/false, ctx, env, a_i8);
+  }
+
+  // Multiple kc ranges: exercises MMAddC accumulation and the once-only
+  // application of the u8 bias correction.
+  TestCase<float>(1, 20000, 8, false, ctx, env, a_i8);
+  TestCase<float>(4, 20000, 64, true, ctx, env, a_i8);
+  TestCase<float>(32, 12345, 64, true, ctx, env, a_i8);
+
+  // BF16 output. The tolerance is loose because `MMAddC` accumulates through
+  // `C`, so with several kc ranges the intermediate sums are rounded to BF16;
+  // the control below shows the existing kernel does the same.
+  TestCase<BF16>(4, 1153, 64, false, ctx, env, a_i8);
+  TestCase<BF16>(32, 20000, 64, false, ctx, env, a_i8);
+  TestCase<BF16>(32, 20000, 64, true, ctx, env, a_i8);
+  ControlBF16OutputError(4, 1153, 64, ctx, env);
+  ControlBF16OutputError(32, 20000, 64, ctx, env);
+
+  // Weights with a large nonzero channel mean, across several kc ranges. This
+  // is the case that a whole-K bias correction gets badly wrong.
+  TestCase<float>(32, 20000, 64, true, ctx, env, a_i8, /*b_mean=*/3.0f);
+  TestCase<BF16>(32, 20000, 64, true, ctx, env, a_i8, /*b_mean=*/3.0f);
+}
+
+// NOLINTNEXTLINE(google-readability-namespace-comments)
+}  // namespace HWY_NAMESPACE
+}  // namespace gcpp
+HWY_AFTER_NAMESPACE();
+
+#if HWY_ONCE
+namespace gcpp {
+size_t g_failures = 0;
+HWY_EXPORT(TestAll);
+}  // namespace gcpp
+
+int main(int /*argc*/, char** /*argv*/) {
+  HWY_DYNAMIC_DISPATCH(gcpp::TestAll)();
+  const size_t failures = gcpp::g_failures;
+  printf("%s (%zu failures)\n", failures == 0 ? "PASS" : "FAIL", failures);
+  return failures == 0 ? 0 : 1;
+}
+#endif  // HWY_ONCE
diff --git a/ops/ops-inl.h b/ops/ops-inl.h
index 8e6b0113..c9e0d00c 100644
--- a/ops/ops-inl.h
+++ b/ops/ops-inl.h
@@ -53,6 +53,7 @@
 
 #include "compression/compress-inl.h"
 #include "ops/dot-inl.h"
+#include "ops/matmul_i8_model-inl.h"
 #include "ops/matmul_static.h"  // includes highway.h
 #include "ops/sum-inl.h"
 #include "hwy/contrib/algo/transform-inl.h"
@@ -72,6 +73,11 @@ MMPerKey* CallMatMul(const MatPtrT<TA>& A, const MatPtr& B,
                      const float* HWY_RESTRICT add, MatMulEnv& env,
                      MatPtrT<TC>& C, const MMOptions& options = MMOptions()) {
   return CallUpcasted(&B, [&](const auto* B_t) {
+    // Experiment: route through the W8A8 kernel if enabled for this tensor.
+    // Returns nullptr when disabled or ineligible, see `matmul_i8_model-inl.h`.
+    if (MMPerKey* per_key = MaybeMatMulI8(A, *B_t, add, env, C, options)) {
+      return per_key;
+    }
     return MatMulStatic(A, *B_t, add, env, C, options);
   });
 }
@@ -79,6 +85,8 @@ MMPerKey* CallMatMul(const MatPtrT<TA>& A, const MatPtr& B,
 static inline void CallTwoMatMul(const MatPtrT<BF16>& A, const MatPtr& B1,
                                  const MatPtr& B2, MatMulEnv& env,
                                  MatPtrT<BF16>& C, const MMOptions& options) {
+  // Experiment, see `matmul_i8_model-inl.h`; nullptr means not enabled here.
+  if (MaybeTwoMatMulI8(A, B1, B2, env, C, options) != nullptr) return;
   return CallUpcastedSame(&B1, &B2, [&](const auto* B1_t, const auto* B2_t) {
     return TwoMatMulStatic(A, *B1_t, *B2_t, env, C, options);
   });
diff --git a/util/zones.cc b/util/zones.cc
index 9cd3475f..3bcfb822 100644
--- a/util/zones.cc
+++ b/util/zones.cc
@@ -177,6 +177,8 @@ const char* CallerName(Callers caller) {
       return "MM.ClusterForN";
     case Callers::kMMClusterForSFC:
       return "MM.ClusterForSFC";
+    case Callers::kMMQuantizeA:
+      return "MM.QuantizeA";
     case Callers::kMMHierForMC:
       return "MM.HierForMC";
     case Callers::kMMHierForMCNC:
diff --git a/util/zones.h b/util/zones.h
index 533182a9..8e426953 100644
--- a/util/zones.h
+++ b/util/zones.h
@@ -101,6 +101,7 @@ enum class Callers {  // Keep sorted
   kMMClusterForMCNC,
   kMMClusterForN,
   kMMClusterForSFC,
+  kMMQuantizeA,
   kMMHierForMC,
   kMMHierForMCNC,
   kMMHierForN,