diff --git a/docs/.nav.yml b/docs/.nav.yml
index eb76761f..9321713a 100644
--- a/docs/.nav.yml
+++ b/docs/.nav.yml
@@ -15,6 +15,7 @@ nav:
       - operators/fused-logp.md
       - operators/linear-logp.md
       - operators/grpo-loss.md
+      - operators/lm_head.md
       - operators/ratio-kl.md
       - operators/sampling.md
       - operators/embedding.md
diff --git a/docs/operators/README.md b/docs/operators/README.md
index a1422758..eb2f5de3 100644
--- a/docs/operators/README.md
+++ b/docs/operators/README.md
@@ -22,6 +22,7 @@ Every operator page should include:
 - [Fused LogP](fused-logp.md)
 - [Fused Linear LogP](linear-logp.md)
 - [GRPO Loss](grpo-loss.md)
+- [LM Head](lm_head.md)
 - [Policy Ratio + KL Penalty](ratio-kl.md)
 - [Sampling](sampling.md)
 - [Token Embedding](embedding.md)
diff --git a/docs/operators/lm_head.md b/docs/operators/lm_head.md
new file mode 100644
index 00000000..a7608385
--- /dev/null
+++ b/docs/operators/lm_head.md
@@ -0,0 +1,127 @@
+# LM Head
+
+The lm_head operator projects hidden states back to vocabulary logits — the final
+layer of the Qwen3/Llama stack. It is a **WS1 ground-truth reference** (issue #108):
+a pure-PyTorch definition of the "correct answer" that downstream fused CUDA/Triton
+kernels are validated against.
+
+- **LM Head** (`NativeLMHeadOp`): `out = hidden @ weight.t() (+ bias)`.
+
+For Qwen3-8B the weight is the output projection `[vocab=151936, hidden=4096]` in the
+HF `nn.Linear` `[out, in]` convention, so it is transposed internally. It is
+**independent** from the embedding table (`tie_word_embeddings=false`) — the two
+weights are not shared — and Qwen3 has **no bias** (`bias=None`).
+
+## Entry Point
+```python
+from rl_engine.kernels.registry import kernel_registry
+
+lm_head = kernel_registry.get_op("lm_head")
+
+logits = lm_head(hidden, weight)          # [B, S, hidden], [vocab, hidden] -> [B, S, vocab]
+logits = lm_head(hidden, weight, bias=b)  # optional [vocab] bias
+```
+
+The op exposes the WS1 dual-path contract:
+
+- `forward(...)` — projects in the input dtype, returns the input dtype (Axis-B accuracy
+  candidate / dtype-behavior path).
+- `forward_fp32(...)` — upcasts to fp32, accumulates in fp32, returns fp32 (the
+  ground-truth golden path). The matmul runs with autocast disabled and CUDA TF32
+  turned off, so it stays a true fp32 reference regardless of the caller's ambient
+  precision context (the global `allow_tf32` flag is saved and restored around it).
+
+## Backends
+
+| Backend | Wrapper | Native symbol | Status |
+| --- | --- | --- | --- |
+| PyTorch fallback | `NativeLMHeadOp` | None | fp32 ground-truth reference; CPU and any GPU. |
+| CUDA / ROCm / Triton | — | — | Planned: downstream fused kernels validate against this reference. |
+
+## Tensor Contract
+
+| Argument | Shape | Dtype | Requirements |
+| --- | --- | --- | --- |
+| `hidden` | `[B, S, hidden]` (any leading dims) | float (fp16/bf16/fp32) | Hidden states (Qwen3-8B `hidden=4096`). |
+| `weight` | `[vocab, hidden]` | float | Output projection in HF `[out, in]` layout; transposed internally. Qwen3-8B `[151936, 4096]`. |
+| `bias` | `[vocab]` or `None` | float | Optional; Qwen3 has none (`None`). |
+| output | `hidden.shape[:-1] + (vocab,)` | `forward`: hidden dtype · `forward_fp32`: float32 | Logits. |
+
+Output dtype follows `hidden`. Pure function — no randomness, no in-place mutation,
+device/dtype follow the inputs.
+
+> **Difference from the bare `matmul` op**: lm_head takes the weight in HF `[out, in]`
+> layout and transposes it internally (`weight.t()`); the `matmul` op computes a bare
+> `a @ b` with no transpose. Do not use them interchangeably.
+
+## Dispatch Behavior
+
+`kernel_registry.get_op("lm_head")` resolves through the `OpBackend` priority map. On
+`cuda` / `rocm` / `cpu` the only registered backend today is the PyTorch native op
+(`PYTORCH_NATIVE_LM_HEAD`), so every device dispatches to the fp32 reference. When fused
+kernels land, they are prepended to the priority list and the native op becomes the fallback.
+
+## Accuracy
+
+Reference semantics (`forward_fp32`):
+
+```python
+out = hidden.float() @ weight.float().t()
+if bias is not None:
+    out = out + bias.float()
+```
+
+- **Ground truth**: `forward_fp32` accumulates in and returns fp32, with autocast and
+  CUDA TF32 disabled so it is a true fp32 reference even if the caller has TF32 or
+  autocast enabled.
+- **Dtype path**: `forward` runs the projection in the input dtype. Because this is a
+  reduction over `hidden`, low-precision accumulation **drifts** from the fp32 reference
+  (unlike the lossless embedding gather). Unlike `forward_fp32`, this path intentionally
+  follows the ambient precision context (it is the dtype-behavior path): with an fp32
+  input it is bitwise-equal to the ground truth **when ambient TF32/autocast is off**, but
+  on a TF32-enabled GPU it tracks real hardware behavior and may drift. bf16/fp16 are
+  always checked with a tolerance.
+- **Axis-B — accuracy tolerance**: measured as max absolute error relative to the output
+  peak magnitude. On a SMALL load point with the real `hidden=4096` reduction length,
+  bf16 drifts ~0.3–0.4% of peak and fp16 ~0.05%. Elementwise `rtol` is not used: many
+  logits are near zero while the accumulated error tracks the reduction length, not the
+  output value.
+- **Axis-A — batch invariance**: a row's logits are independent of the rest of the batch,
+  so the output is bitwise-identical regardless of batch size or padding (`torch.equal`,
+  `atol=0`) — **provided the reduction order is fixed**. Multi-threaded CPU GEMM splits
+  the `hidden` reduction across threads by the `M = batch*seq` dimension, which silently
+  breaks bitwise batch invariance for large `hidden`; the tests pin a single thread to fix
+  the order. On GPU, cuBLAS likewise splits K by `M`, so a bitwise batch-invariant GEMM is
+  a downstream kernel concern, not a free property of `torch.matmul`.
+
+## Performance Notes
+
+Reference operator — no fused kernel or benchmark yet. Downstream fused kernels carry their
+own benchmarks and are measured against this reference for correctness.
+
+## Tests
+
+```bash
+python -m pytest tests/test_lm_head.py -v
+```
+
+Covers: fp32 correctness vs naive matmul (bitwise, with ambient TF32 pinned off),
+`forward_fp32` precision-context safety (true fp32 under ambient autocast + restores the
+global TF32 flag on CPU; numerically beats a TF32 matmul on GPU), bf16/fp16 dtype-path
+accuracy (relative-to-peak tolerance, with `bias`), output shape, bias semantics, Axis-A
+batch invariance (slice + padding, single-thread reduction, all dtypes), input purity,
+gradient flow to `hidden`/`weight` (closed-form check), registry dispatch, and a GPU-only
+smoke test at the real Qwen3-8B dims (`vocab=151936, hidden=4096`) that skips when CUDA or
+GPU memory is unavailable.
+
+## Implementation Files
+
+- `rl_engine/kernels/ops/pytorch/linear/lm_head.py`
+- `rl_engine/kernels/registry.py`
+- `tests/test_lm_head.py`
+
+## Known Limitations
+
+- PyTorch fallback only; no fused CUDA/Triton backend yet (downstream work).
+- Axis-A bitwise batch invariance holds only with a fixed reduction order (single-thread on
+  CPU); a batch-invariant GEMM on GPU is a downstream concern.
diff --git a/rl_engine/kernels/ops/pytorch/linear/lm_head.py b/rl_engine/kernels/ops/pytorch/linear/lm_head.py
new file mode 100644
index 00000000..4604fee7
--- /dev/null
+++ b/rl_engine/kernels/ops/pytorch/linear/lm_head.py
@@ -0,0 +1,124 @@
+# SPDX-License-Identifier: Apache-2.0
+# Copyright (c) 2026 RL-Kernel Contributors
+
+from __future__ import annotations
+
+from contextlib import contextmanager
+from typing import Optional
+
+import torch
+
+
+class NativeLMHeadOp:
+    """
+    Pure PyTorch native language-model-head reference.
+    out = hidden @ weight.t() (+ bias)
+
+    Projects hidden states back to vocabulary logits -- the final layer of the
+    Qwen3/Llama stack. For Qwen3-8B the weight is the output projection
+    ``[vocab=151936, hidden=4096]`` in the HF ``nn.Linear`` ``[out, in]``
+    convention, so it is transposed internally (``weight.t()``). This is the
+    one difference from the bare ``matmul`` op (which computes ``a @ b`` with no
+    transpose) -- do not use them interchangeably. The lm_head weight is
+    *independent* from the embedding table (``tie_word_embeddings=false``), and
+    Qwen3 has no bias (pass ``bias=None``).
+
+    Unlike embedding (a lossless row gather), this is a reduction over the
+    ``hidden`` dimension: the low-precision ``forward`` path accumulates in the
+    input dtype and therefore drifts from the fp32 ``forward_fp32`` ground
+    truth. Axis-B accuracy uses a tolerance (``torch.allclose``), not bitwise
+    equality. Axis-A batch invariance still holds bitwise within a single dtype
+    (each output row reduces over ``hidden`` independently of the batch).
+    """
+
+    def __init__(self) -> None:
+        """No state; the op is a pure function over (hidden, weight, bias)."""
+
+    def __call__(
+        self,
+        hidden: torch.Tensor,
+        weight: torch.Tensor,
+        *,
+        bias: Optional[torch.Tensor] = None,
+    ) -> torch.Tensor:
+        """Alias for ``forward`` so the op is callable like a module."""
+        return self.forward(hidden, weight, bias=bias)
+
+    def forward(
+        self,
+        hidden: torch.Tensor,
+        weight: torch.Tensor,
+        *,
+        bias: Optional[torch.Tensor] = None,
+    ) -> torch.Tensor:
+        """
+        Canonical entry: project in the input dtype, output the input dtype.
+        This is the dtype-behavior path used as the Axis-B accuracy candidate.
+        """
+        return self._lm_head(
+            hidden, weight, bias, compute_dtype=hidden.dtype, output_dtype=hidden.dtype
+        )
+
+    def forward_fp32(
+        self,
+        hidden: torch.Tensor,
+        weight: torch.Tensor,
+        *,
+        bias: Optional[torch.Tensor] = None,
+    ) -> torch.Tensor:
+        """Ground truth: upcast to fp32, accumulate in fp32, force fp32 output.
+
+        The matmul is wrapped to disable autocast and TF32 so this stays a true
+        fp32 reference regardless of the caller's ambient precision context.
+        """
+        return self._lm_head(
+            hidden,
+            weight,
+            bias,
+            compute_dtype=torch.float32,
+            output_dtype=torch.float32,
+            strict_fp32=True,
+        )
+
+    # ------------------------------------------------------------------ #
+    # Helpers
+    # ------------------------------------------------------------------ #
+    @staticmethod
+    def _lm_head(
+        hidden: torch.Tensor,
+        weight: torch.Tensor,
+        bias: Optional[torch.Tensor],
+        *,
+        compute_dtype: torch.dtype,
+        output_dtype: torch.dtype,
+        strict_fp32: bool = False,
+    ) -> torch.Tensor:
+        """Core matmul: cast to ``compute_dtype``, project, optionally add bias, cast out.
+
+        When ``strict_fp32`` is set, the matmul runs with autocast disabled and
+        CUDA TF32 turned off so the fp32 reference path is not silently downcast
+        by the caller's ambient autocast/TF32 settings.
+        """
+        h = hidden.to(compute_dtype)
+        w = weight.to(compute_dtype)
+        # [..., hidden] @ [hidden, vocab] -> [..., vocab]; weight is [vocab, hidden] (HF [out, in]).
+        if strict_fp32:
+            with NativeLMHeadOp._strict_fp32_matmul(h.device.type):
+                out = h @ w.t()
+        else:
+            out = h @ w.t()
+        if bias is not None:
+            out = out + bias.to(compute_dtype)
+        return out.to(output_dtype)
+
+    @staticmethod
+    @contextmanager
+    def _strict_fp32_matmul(device_type: str):
+        """Disable autocast and TF32 for a true fp32 matmul, restoring state after."""
+        prev_tf32 = torch.backends.cuda.matmul.allow_tf32
+        torch.backends.cuda.matmul.allow_tf32 = False
+        try:
+            with torch.autocast(device_type=device_type, enabled=False):
+                yield
+        finally:
+            torch.backends.cuda.matmul.allow_tf32 = prev_tf32
diff --git a/rl_engine/kernels/registry.py b/rl_engine/kernels/registry.py
index d658ab20..5d866403 100644
--- a/rl_engine/kernels/registry.py
+++ b/rl_engine/kernels/registry.py
@@ -56,6 +56,8 @@ class OpBackend(Enum, metaclass=_KernelEnumMeta):
     PYTORCH_NATIVE_SILU = "rl_engine.kernels.ops.pytorch.activation.swiglu.NativeSiLUOp"
     PYTORCH_NATIVE_SWIGLU = "rl_engine.kernels.ops.pytorch.activation.swiglu.NativeSwiGLUOp"
 
+    # WS1 pure-PyTorch ground-truth linear ops
+    PYTORCH_NATIVE_LM_HEAD = "rl_engine.kernels.ops.pytorch.linear.lm_head.NativeLMHeadOp"
     # WS1 pure-PyTorch ground-truth embedding ops
     PYTORCH_NATIVE_EMBEDDING = "rl_engine.kernels.ops.pytorch.linear.embedding.NativeEmbeddingOp"
 
@@ -94,6 +96,7 @@ def __init__(self):
                 "grpo_loss": [OpBackend.TRITON_GRPO_LOSS, OpBackend.PYTORCH_GRPO_LOSS],
                 "linear_logp": [OpBackend.TRITON_LINEAR_LOGP, OpBackend.PYTORCH_LINEAR_LOGP],
                 "ratio_kl": [OpBackend.TRITON_RATIO_KL, OpBackend.PYTORCH_RATIO_KL],
+                "lm_head": [OpBackend.PYTORCH_NATIVE_LM_HEAD],
                 "embedding": [OpBackend.PYTORCH_NATIVE_EMBEDDING],
                 "silu": [OpBackend.PYTORCH_NATIVE_SILU],
                 "swiglu": [OpBackend.PYTORCH_NATIVE_SWIGLU],
@@ -109,6 +112,7 @@ def __init__(self):
                 "grpo_loss": [OpBackend.TRITON_GRPO_LOSS, OpBackend.PYTORCH_GRPO_LOSS],
                 "linear_logp": [OpBackend.TRITON_LINEAR_LOGP, OpBackend.PYTORCH_LINEAR_LOGP],
                 "ratio_kl": [OpBackend.TRITON_RATIO_KL, OpBackend.PYTORCH_RATIO_KL],
+                "lm_head": [OpBackend.PYTORCH_NATIVE_LM_HEAD],
                 "embedding": [OpBackend.PYTORCH_NATIVE_EMBEDDING],
                 "silu": [OpBackend.PYTORCH_NATIVE_SILU],
                 "swiglu": [OpBackend.PYTORCH_NATIVE_SWIGLU],
@@ -119,6 +123,7 @@ def __init__(self):
                 "grpo_loss": [OpBackend.PYTORCH_GRPO_LOSS],
                 "linear_logp": [OpBackend.PYTORCH_LINEAR_LOGP],
                 "ratio_kl": [OpBackend.PYTORCH_RATIO_KL],
+                "lm_head": [OpBackend.PYTORCH_NATIVE_LM_HEAD],
                 "embedding": [OpBackend.PYTORCH_NATIVE_EMBEDDING],
                 "silu": [OpBackend.PYTORCH_NATIVE_SILU],
                 "swiglu": [OpBackend.PYTORCH_NATIVE_SWIGLU],
diff --git a/tests/test_lm_head.py b/tests/test_lm_head.py
new file mode 100644
index 00000000..15571e2b
--- /dev/null
+++ b/tests/test_lm_head.py
@@ -0,0 +1,337 @@
+# SPDX-License-Identifier: Apache-2.0
+# Copyright (c) 2026 RL-Kernel Contributors
+"""Tests for NativeLMHeadOp (ISSUE #108 WS1 ground-truth baseline).
+
+The lm_head projects hidden states to vocab logits: out = hidden @ weight.t()
+(+ bias). Unlike embedding (a lossless gather), this is a *reduction* over the
+hidden dimension, so:
+
+  * Axis-B (accuracy): the low-precision ``forward`` path accumulates in the
+    input dtype and drifts from the fp32 ``forward_fp32`` ground truth. It is
+    checked with a tolerance (relative to the output peak magnitude), not
+    bitwise -- elementwise rtol is useless here because many logits are near
+    zero while the accumulated error tracks the reduction length, not the
+    output value.
+  * Axis-A (batch invariance): still bitwise within a single dtype, but only
+    once the CPU reduction order is pinned. Multi-threaded CPU GEMM splits the
+    K (=hidden) reduction across threads differently depending on the M
+    (=batch*seq) dimension, which silently breaks bitwise batch invariance for
+    large hidden. ``_single_thread`` fixes the reduction order; this is the
+    local stand-in for the planned testing/determinism.py::deterministic_context.
+"""
+
+import contextlib
+
+import pytest
+import torch
+
+from rl_engine.kernels.ops.pytorch.linear.lm_head import NativeLMHeadOp
+from rl_engine.kernels.registry import kernel_registry
+
+# Qwen3-8B architecture (synthetic tensors, no weight download). The vocab is
+# shrunk -- it is just the number of independent output dot products, so the
+# logic is identical at any vocab. The reduction dim (hidden) is kept at the
+# *real* 4096 so the Axis-B drift is representative: weight is only
+# [vocab=128, hidden=4096] ~ 2 MB, trivial on CPU. The full vocab is exercised
+# by the GPU smoke test below.
+_VOCAB = 128  # shrunk; real value: _QWEN3_VOCAB
+_HIDDEN = 4096  # real Qwen3-8B reduction dim -- kept real for representative drift
+
+# Real Qwen3-8B output-projection dims: 151936 x 4096 ~ 2.49 GB in fp32.
+_QWEN3_VOCAB = 151936
+_QWEN3_HIDDEN = 4096
+
+# Axis-B: max abs error as a fraction of the output peak magnitude. Calibrated
+# from measured SMALL drift (bf16 ~0.3% of peak, fp16 ~0.04%) with headroom.
+_DTYPE_REL_PEAK = {torch.bfloat16: 1.0e-2, torch.float16: 2.0e-3}
+
+
+def _cpu_fp16_matmul_supported() -> bool:
+    """Probe whether this CPU backend implements float16 matmul."""
+    try:
+        _ = torch.randn(2, 2, dtype=torch.float16) @ torch.randn(2, 2, dtype=torch.float16)
+        return True
+    except RuntimeError:
+        return False
+
+
+# CPU half-precision matmul is backend/ISA-dependent (AVX512_FP16, AMX) and may
+# be unimplemented on some runners -- gate the fp16 axis so a missing kernel
+# skips rather than fails the test.
+_FP16_IF_CPU_MATMUL_SUPPORTED = pytest.param(
+    torch.float16,
+    marks=pytest.mark.skipif(
+        not _cpu_fp16_matmul_supported(),
+        reason="CPU float16 matmul unsupported on this backend",
+    ),
+)
+_DTYPES_AXIS_B = (torch.bfloat16, _FP16_IF_CPU_MATMUL_SUPPORTED)
+_DTYPES_AXIS_A = (torch.float32, torch.bfloat16, _FP16_IF_CPU_MATMUL_SUPPORTED)
+
+
+@contextlib.contextmanager
+def _single_thread():
+    """Pin CPU GEMM to one thread so the K reduction order is M-independent."""
+    prev = torch.get_num_threads()
+    torch.set_num_threads(1)
+    try:
+        yield
+    finally:
+        torch.set_num_threads(prev)
+
+
+# Shared helpers -- fixed-seed Generator for determinism / reproducibility.
+def _rand_hidden(batch, seq, hidden=_HIDDEN, *, seed, dtype=torch.float32):
+    """Fixed-seed random hidden-state tensor for reproducibility."""
+    gen = torch.Generator().manual_seed(seed)
+    return torch.randn(batch, seq, hidden, generator=gen, dtype=dtype)
+
+
+def _rand_weight(vocab=_VOCAB, hidden=_HIDDEN, *, seed, dtype=torch.float32):
+    """Fixed-seed random lm_head weight tensor for reproducibility."""
+    gen = torch.Generator().manual_seed(seed)
+    return torch.randn(vocab, hidden, generator=gen, dtype=dtype)
+
+
+# Correctness of the fp32 ground truth: forward_fp32 == naive fp32 matmul,
+# bitwise. forward_fp32 disables autocast/TF32, so it is a true fp32 reference
+# unconditionally. The fp32 dtype path (forward) follows the ambient precision
+# context, so it is bitwise-equal to the ground truth only when TF32/autocast is
+# off (the default here); only bf16/fp16 introduce drift.
+def test_native_lm_head_fp32_matches_naive_matmul():
+    """forward_fp32 (and fp32 forward, TF32 off) is bitwise-equal to a naive fp32 matmul."""
+    hidden = _rand_hidden(2, 5, seed=1)
+    weight = _rand_weight(seed=1)
+
+    # Pin TF32 off so the fp32 forward path and the naive reference below are both
+    # true fp32 regardless of the machine's global default, making this assertion
+    # deterministic. forward_fp32 disables TF32 internally and does not need this.
+    prev_tf32 = torch.backends.cuda.matmul.allow_tf32
+    torch.backends.cuda.matmul.allow_tf32 = False
+    try:
+        naive = hidden.float() @ weight.float().t()
+        assert torch.equal(NativeLMHeadOp().forward_fp32(hidden, weight), naive)
+        # fp32 forward path computes in fp32 too -> bitwise equal to ground truth.
+        assert torch.equal(NativeLMHeadOp().forward(hidden, weight), naive)
+    finally:
+        torch.backends.cuda.matmul.allow_tf32 = prev_tf32
+
+
+def test_forward_fp32_ignores_ambient_autocast_and_restores_tf32():
+    """forward_fp32 is a strict fp32 reference under ambient autocast/TF32 settings."""
+    op = NativeLMHeadOp()
+    hidden = _rand_hidden(2, 5, hidden=32, seed=11)
+    weight = _rand_weight(vocab=7, hidden=32, seed=11)
+    gen = torch.Generator().manual_seed(11)
+    bias = torch.randn(7, generator=gen)
+    ref = hidden.float() @ weight.float().t() + bias.float()
+
+    prev_tf32 = torch.backends.cuda.matmul.allow_tf32
+    torch.backends.cuda.matmul.allow_tf32 = True
+    try:
+        with torch.autocast(device_type="cpu", dtype=torch.bfloat16):
+            out = op.forward_fp32(hidden, weight, bias=bias)
+
+        assert out.dtype == torch.float32
+        assert torch.equal(out, ref)
+        assert torch.backends.cuda.matmul.allow_tf32 is True
+    finally:
+        torch.backends.cuda.matmul.allow_tf32 = prev_tf32
+
+
+# GPU counterpart of the test above: TF32 only changes the numbers on CUDA, so
+# the "TF32 is disabled" half of forward_fp32 can only be exercised on a GPU.
+# With TF32 forced on globally, forward_fp32 must stay true fp32 (tight against a
+# high-precision fp64 reference) and never be worse than a plain TF32 matmul. The
+# assertion holds whether or not the GPU actually uses TF32 for fp32 matmul.
+@pytest.mark.skipif(not torch.cuda.is_available(), reason="needs a CUDA GPU to exercise TF32")
+def test_forward_fp32_disables_tf32_on_gpu():
+    """On a TF32-enabled GPU, forward_fp32 stays true fp32, no worse than a TF32 matmul."""
+    device = torch.device("cuda")
+    gen = torch.Generator(device=device).manual_seed(21)
+    hidden = torch.randn(2, 16, _HIDDEN, generator=gen, dtype=torch.float32, device=device)
+    weight = torch.randn(256, _HIDDEN, generator=gen, dtype=torch.float32, device=device)
+    ref = (hidden.double() @ weight.double().t()).float()  # high-precision reference
+
+    prev_tf32 = torch.backends.cuda.matmul.allow_tf32
+    torch.backends.cuda.matmul.allow_tf32 = True  # hostile ambient setting
+    try:
+        strict = NativeLMHeadOp().forward_fp32(hidden, weight)  # must ignore TF32
+        tf32 = hidden @ weight.t()  # plain matmul -> TF32 if the GPU supports it
+    finally:
+        torch.backends.cuda.matmul.allow_tf32 = prev_tf32
+
+    peak = ref.abs().max().item()
+    strict_err = (strict - ref).abs().max().item()
+    tf32_err = (tf32 - ref).abs().max().item()
+    print(
+        f"\n[lm_head fp32-vs-tf32] strict_err={strict_err:.3g} "
+        f"tf32_err={tf32_err:.3g} peak={peak:.3g}"
+    )
+    # forward_fp32 is fp32-tight (well under the TF32 ~10-bit-mantissa drift floor)...
+    assert strict_err <= 1.0e-3 * peak
+    # ...and never worse than a plain matmul (far tighter when the GPU uses TF32).
+    assert strict_err <= tf32_err
+
+
+# Axis-B accuracy: the low-precision dtype path drifts from the fp32 reference
+# by a bounded fraction of the output peak. Errors/stats are printed for the PR.
+@pytest.mark.parametrize("dtype", _DTYPES_AXIS_B)
+def test_native_lm_head_dtype_path_accuracy(dtype: torch.dtype):
+    """Axis-B: the low-precision path drifts from fp32 by a bounded fraction of the output peak."""
+    op = NativeLMHeadOp()
+    hidden = _rand_hidden(2, 16, seed=2)
+    weight = _rand_weight(seed=2)
+
+    ref = op.forward_fp32(hidden, weight)  # fp32 ground truth
+    cand = op.forward(hidden.to(dtype), weight.to(dtype))  # dtype path
+    assert cand.dtype == dtype
+
+    err = (cand.float() - ref).abs()
+    peak = ref.abs().max()
+    max_abs, mean_abs = err.max().item(), err.mean().item()
+    print(f"\n[lm_head {dtype}] max_abs={max_abs:.4g} mean_abs={mean_abs:.4g} peak={peak:.4g}")
+    assert max_abs <= _DTYPE_REL_PEAK[dtype] * peak.item()
+
+
+# Output shape must be hidden.shape[:-1] + (vocab,).
+def test_native_lm_head_output_shape():
+    """Output shape is hidden.shape[:-1] + (vocab,)."""
+    hidden = _rand_hidden(3, 7, seed=3)
+    weight = _rand_weight(seed=3)
+    out = NativeLMHeadOp().forward(hidden, weight)
+    assert out.shape == (3, 7, _VOCAB)
+
+
+# Bias: None (Qwen3 default) is a plain matmul; a provided [vocab] bias is added.
+def test_native_lm_head_bias():
+    """bias=None is a plain matmul; a [vocab] bias is added elementwise."""
+    op = NativeLMHeadOp()
+    hidden = _rand_hidden(2, 4, seed=4)
+    weight = _rand_weight(seed=4)
+    gen = torch.Generator().manual_seed(4)
+    bias = torch.randn(_VOCAB, generator=gen)
+
+    no_bias = op.forward_fp32(hidden, weight)
+    with_bias = op.forward_fp32(hidden, weight, bias=bias)
+    assert torch.equal(with_bias, no_bias + bias.float())
+    # default is bias=None (== no bias term).
+    assert torch.equal(op.forward_fp32(hidden, weight, bias=None), no_bias)
+
+
+# Axis A -- batch invariance, bitwise (the WS1 "aligned" property). A row's
+# logits must not depend on how many other rows share the batch. Compute on the
+# full input once, then slice -- never compute a slice on its own. Requires the
+# pinned single-thread reduction order (see module docstring).
+@pytest.mark.parametrize("dtype", _DTYPES_AXIS_A)
+def test_lm_head_batch_invariance_slice(dtype: torch.dtype):
+    """Axis-A: a row's logits are bitwise-independent of how many rows share the batch."""
+    op = NativeLMHeadOp()
+    hidden = _rand_hidden(8, 32, seed=5).to(dtype)
+    weight = _rand_weight(seed=5).to(dtype)
+    with _single_thread():
+        full = op.forward(hidden, weight)  # compute on full batch...
+        assert torch.equal(op.forward(hidden[:1], weight), full[:1])  # ...then slice
+        assert torch.equal(op.forward(hidden[3:5], weight), full[3:5])
+
+
+@pytest.mark.parametrize("dtype", _DTYPES_AXIS_A)
+def test_lm_head_batch_invariance_with_padding(dtype: torch.dtype):
+    """Padding extra seq positions must not perturb the real ones (bitwise)."""
+    op = NativeLMHeadOp()
+    weight = _rand_weight(seed=6).to(dtype)
+    real = _rand_hidden(4, 10, seed=6).to(dtype)
+    pad = _rand_hidden(4, 6, seed=99).to(dtype)  # 6 extra padding positions
+    padded = torch.cat([real, pad], dim=1)  # concat along seq
+    with _single_thread():
+        assert torch.equal(op.forward(padded, weight)[:, :10], op.forward(real, weight))
+
+
+# Purity -- no input may be mutated in place.
+def test_lm_head_inputs_not_mutated():
+    """Purity: no input tensor is mutated in place."""
+    op = NativeLMHeadOp()
+    hidden = _rand_hidden(2, 8, seed=7)
+    weight = _rand_weight(seed=7)
+    gen = torch.Generator().manual_seed(7)
+    bias = torch.randn(_VOCAB, generator=gen)
+    h_c, w_c, b_c = hidden.clone(), weight.clone(), bias.clone()
+    op.forward(hidden, weight, bias=bias)
+    op.forward_fp32(hidden, weight, bias=bias)
+    assert torch.equal(hidden, h_c) and torch.equal(weight, w_c) and torch.equal(bias, b_c)
+
+
+# Gradient (fp32 autograd = backward golden source). Backprop a *random*
+# cotangent, not out.sum(): an all-ones cotangent collapses both grads to column
+# sums (every row of each grad identical), so a transposed / mis-contracted
+# backward could still pass under that symmetry. A random dy exercises the full
+# contraction. For out = hidden @ weight.t() (weight is HF [V, K]) the exact
+# closed forms are dL/dhidden = dy @ weight and dL/dweight = dy^T @ hidden.
+def test_lm_head_gradient_flows():
+    """fp32 autograd matches the closed-form grads under a random cotangent."""
+    op = NativeLMHeadOp()
+    hidden = _rand_hidden(2, 4, seed=8).requires_grad_(True)
+    weight = _rand_weight(seed=8).requires_grad_(True)
+    out = op.forward_fp32(hidden, weight)
+
+    gen = torch.Generator().manual_seed(8)
+    dy = torch.randn(out.shape, generator=gen, dtype=out.dtype)
+    out.backward(dy)
+
+    assert torch.isfinite(hidden.grad).all() and torch.isfinite(weight.grad).all()
+    assert hidden.grad.shape == hidden.shape and weight.grad.shape == weight.shape
+    exp_h = dy @ weight.detach()  # [.., V] @ [V, K] -> [.., K]
+    exp_w = dy.reshape(-1, _VOCAB).t() @ hidden.detach().reshape(-1, _HIDDEN)  # [V, K]
+    torch.testing.assert_close(hidden.grad, exp_h, rtol=1e-5, atol=1e-5)
+    torch.testing.assert_close(weight.grad, exp_w, rtol=1e-5, atol=1e-5)
+
+
+# Registry dispatch -- "lm_head" resolves to NativeLMHeadOp.
+def test_registry_dispatches_native_lm_head_op():
+    """The registry resolves "lm_head" to NativeLMHeadOp."""
+    assert isinstance(kernel_registry.get_op("lm_head"), NativeLMHeadOp)
+
+
+# --------------------------------------------------------------------------- #
+# Qwen3-8B real-shape smoke test
+# --------------------------------------------------------------------------- #
+# Exercises the real output-projection dims (vocab=151936, hidden=4096). The
+# fp32 weight is ~2.5 GB, so this is GPU-only and skips when CUDA is absent or
+# there is not enough free memory. The shrunk-vocab tests above already cover
+# the logic; this validates the real vocab width and hidden reduction length at
+# a small (batch, seq) load point.
+def _enough_gpu_memory(num_bytes: int) -> bool:
+    """Return True only if CUDA is present and has free memory with headroom."""
+    if not torch.cuda.is_available():
+        return False
+    try:
+        free, _ = torch.cuda.mem_get_info()
+    except RuntimeError:
+        return False
+    return free > int(num_bytes * 1.5)  # headroom for the logits output
+
+
+@pytest.mark.skipif(
+    not _enough_gpu_memory(_QWEN3_VOCAB * _QWEN3_HIDDEN * 4),
+    reason="needs a CUDA GPU with room for the ~2.5 GB fp32 Qwen3-8B lm_head weight",
+)
+def test_native_lm_head_qwen3_8b_real_shape():
+    """GPU smoke test at real Qwen3-8B dims (vocab=151936, hidden=4096)."""
+    device = torch.device("cuda")
+    op = NativeLMHeadOp()
+
+    # SMALL load point (batch=2, seq=16) at the real model dims.
+    gen = torch.Generator(device=device).manual_seed(0)
+    hidden = torch.randn(2, 16, _QWEN3_HIDDEN, generator=gen, dtype=torch.float32, device=device)
+    weight = torch.randn(
+        _QWEN3_VOCAB, _QWEN3_HIDDEN, generator=gen, dtype=torch.float32, device=device
+    )
+
+    out = op.forward_fp32(hidden, weight)
+    assert out.shape == (2, 16, _QWEN3_VOCAB)
+    assert out.dtype == torch.float32
+    # Bitwise equal to the naive fp32 matmul (same call, same inputs).
+    assert torch.equal(out, hidden @ weight.t())
+    # NB: Axis-A bitwise is asserted on CPU (single-thread reduction) above.
+    # On GPU it is NOT free -- cuBLAS also splits K by the M dimension, so a
+    # batch-invariant GEMM is a downstream kernel concern, not validated here.