diff --git a/python/freetoken/kernel/gluon/nvfp4/moe/sm120_w4a4.py b/python/freetoken/kernel/gluon/nvfp4/moe/sm120_w4a4.py
new file mode 100644
index 000000000..2a2c1430c
--- /dev/null
+++ b/python/freetoken/kernel/gluon/nvfp4/moe/sm120_w4a4.py
@@ -0,0 +1,536 @@
+"""NVFP4 W4A4 MoE prefill GEMMs for sm_120: NVFP4 weights times NVFP4 activations on mma.sync kind::mxf4nvf4.
+
+The activations are quantized under static global scales (the checkpoint's input_scale), so each GEMM input is e2m1 codes, per-16 e4m3 scales and one fp32 scalar.
+:func:`gate_up_gemm` applies the gated activation in its epilogue and writes the down GEMM's NVFP4 input; :func:`down_gemm` scales by the routing weight and adds each route into its token's output row with bf16 atomics (the sum order, and so the result, varies from run to run).
+Routing follows the W4A16 kernels (moe_align sorted ids, expert per M block, swap-AB accumulator [BN, BM]); :mod:`freetoken.kernel.nvfp4_moe` picks the constexprs.
+"""
+
+import torch
+import triton
+from triton.experimental import gluon
+from triton.experimental.gluon import language as gl
+from triton.experimental.gluon.language.nvidia.ampere import async_copy
+from triton.experimental.gluon.language.nvidia.hopper import mbarrier, tma
+from triton.experimental.gluon.nvidia.hopper import TensorDescriptor
+
+from ..utils import mma_nvf4, scale_layout
+
+# gated activations of the gate_up epilogue, as gated_act_and_mul computes them
+ACTIVATIONS = {"silu": 0, "gelu_tanh": 1, "swigluoai": 2, "swiglu_clamp": 3}
+
+SCALE_SMEM_LAYOUT = gl.constexpr(gl.SwizzledSharedLayout(1, 1, 1, [1, 0]))
+
+
+@gluon.jit
+def _tile(t, num_pid_m, num_pid_n, GROUP_SIZE_M: gl.constexpr, N_GROUP: gl.constexpr):
+    if N_GROUP > 0:
+        # N-outer strips of N_GROUP tiles: concurrent CTAs add into a narrow output strip that stays in L2
+        per_strip = num_pid_m * N_GROUP
+        strip = t // per_strip
+        rem = t % per_strip
+        width = min(N_GROUP, num_pid_n - strip * N_GROUP)
+        return rem // width, strip * N_GROUP + rem % width
+    num_pid_in_group = GROUP_SIZE_M * num_pid_n
+    first_pid_m = (t // num_pid_in_group) * GROUP_SIZE_M
+    group_size_m = min(num_pid_m - first_pid_m, GROUP_SIZE_M)
+    return first_pid_m + ((t % num_pid_in_group) % group_size_m), (t % num_pid_in_group) // group_size_m
+
+
+@gluon.jit
+def _scale_ptrs(ws_ptr, xs_ptr, ids_ptr, w_row, w_row_up, x_div, pid_m, num_valid, ks4, n_rows,
+                BM: gl.constexpr, BN: gl.constexpr, K_TAIL: gl.constexpr, GATED: gl.constexpr, NWL: gl.constexpr):
+    # scale rows are K / 16 bytes: 16-byte copies when K % 256 == 0, else 4-byte copies zeroed past the row end
+    SV: gl.constexpr = 4 if K_TAIL else 16
+    sll: gl.constexpr = gl.BlockedLayout([1, SV], [32 * SV // 16, 16 // SV], [NWL, 1], [1, 0])
+    s_cols = gl.arange(0, 16, layout=gl.SliceLayout(0, sll))
+    # both forms equal K // 16; the second tells the compiler the row is a multiple of 16 bytes
+    row_len = ks4 * 4 if K_TAIL else (ks4 // 4) * 16
+    r = gl.arange(0, BN, layout=gl.SliceLayout(1, sll))
+    if GATED:
+        sw_rows = gl.where(r < BN // 2, w_row + r, w_row_up + r - BN // 2)
+    else:
+        sw_rows = w_row + r
+    sw_ptrs = ws_ptr + sw_rows.to(gl.int64)[:, None] * row_len + s_cols[None, :]
+    sw_ok = (sw_rows < n_rows)[:, None]
+    sx_tok = gl.load(ids_ptr + pid_m * BM + gl.arange(0, BM, layout=gl.SliceLayout(1, sll)))
+    # the x scale rows are contiguous K / 16 bytes too; row_len keeps the divisibility a stride argument loses
+    sx_ptrs = xs_ptr + ((sx_tok // x_div).to(gl.int64) * row_len)[:, None] + s_cols[None, :]
+    sx_ok = (sx_tok < num_valid)[:, None]
+    return sw_ptrs, sw_ok, sx_ptrs, sx_ok, s_cols, row_len
+
+
+@gluon.jit
+def _load_tile(xq_ptr, xs_ptr, ids_ptr, eid_ptr, w_desc, ws_ptr, bufs, bars, pid_m, pid_n, N, nk, it, jt, num_valid,
+               stride_xq, ks4, n_rows, x_div, K_TAIL: gl.constexpr, GATED: gl.constexpr):
+    x_bufs, w_bufs, sw_bufs, sx_bufs = bufs
+    ready, empty, s_ready, s_empty = bars
+    S: gl.constexpr = x_bufs.type.shape[0]
+    BM: gl.constexpr = x_bufs.type.shape[1]
+    KB: gl.constexpr = x_bufs.type.shape[2]
+    BN: gl.constexpr = w_bufs.type.shape[1]
+    SPC: gl.constexpr = 128 // KB
+    NWL: gl.constexpr = gl.num_warps()
+
+    if GATED:
+        # rows [0, BN/2) of the tile are gate rows, [BN/2, BN) the matching up rows N/2 further down
+        w_row = gl.load(eid_ptr + pid_m) * N + pid_n * (BN // 2)
+        w_row_up = w_row + N // 2
+    else:
+        w_row = gl.load(eid_ptr + pid_m) * N + pid_n * BN
+        w_row_up = w_row
+
+    # x rows: 16 bytes per cp.async, KB // 16 lanes per row
+    xl: gl.constexpr = gl.BlockedLayout([1, 16], [32 * 16 // KB, KB // 16], [NWL, 1], [1, 0])
+    tok = gl.load(ids_ptr + pid_m * BM + gl.arange(0, BM, layout=gl.SliceLayout(1, xl)))
+    x_ptrs = xq_ptr + ((tok // x_div).to(gl.int64) * stride_xq)[:, None] + gl.arange(0, KB, layout=gl.SliceLayout(0, xl))[None, :]
+    x_mask = (tok < num_valid)[:, None]
+    sw_ptrs, sw_ok, sx_ptrs, sx_ok, s_cols, row_len = _scale_ptrs(
+        ws_ptr, xs_ptr, ids_ptr, w_row, w_row_up, x_div, pid_m, num_valid, ks4, n_rows, BM, BN, K_TAIL, GATED, NWL)
+
+    for kk in range(nk):
+        if kk % SPC == 0:
+            jl = kk // SPC
+            j = jt + jl
+            ss = j % 2
+            mbarrier.wait(s_empty.index(ss), (j // 2 & 1) ^ 1)
+            col_ok = (jl * 16 + s_cols < row_len)[None, :]
+            async_copy.async_load(sw_bufs.index(ss), sw_ptrs + jl * 16, mask=sw_ok & col_ok)
+            async_copy.async_load(sx_bufs.index(ss), sx_ptrs + jl * 16, mask=sx_ok & col_ok)
+            async_copy.mbarrier_arrive(s_ready.index(ss))
+            mbarrier.arrive(s_ready.index(ss))
+        c = it + kk
+        st = c % S
+        mbarrier.wait(empty.index(st), (c // S & 1) ^ 1)
+        # the cp.async arrive goes first: it raises the pending count before the TMA expect can complete the phase
+        async_copy.async_load(x_bufs.index(st), x_ptrs + kk * KB, mask=x_mask)
+        async_copy.mbarrier_arrive(ready.index(st))
+        if GATED:
+            mbarrier.expect(ready.index(st), 2 * w_desc.block_type.nbytes)
+            tma.async_load(w_desc, [w_row, kk * KB], ready.index(st), w_bufs.index(st).slice(0, BN // 2))
+            tma.async_load(w_desc, [w_row_up, kk * KB], ready.index(st), w_bufs.index(st).slice(BN // 2, BN // 2))
+        else:
+            mbarrier.expect(ready.index(st), w_desc.block_type.nbytes)
+            tma.async_load(w_desc, [w_row, kk * KB], ready.index(st), w_bufs.index(st))
+
+
+@gluon.jit
+def _load_partition(xq_ptr, xs_ptr, ids_ptr, eid_ptr, w_desc, ws_ptr, bufs, bars, N, nk, num_pid_m, num_valid,
+                    stride_xq, ks4, n_rows, x_div, K_TAIL: gl.constexpr, GROUP_SIZE_M: gl.constexpr,
+                    GATED: gl.constexpr, N_GROUP: gl.constexpr):
+    BN: gl.constexpr = bufs[1].type.shape[1]
+    pid_m, pid_n = _tile(gl.program_id(0), num_pid_m, gl.cdiv(N, BN), GROUP_SIZE_M, N_GROUP)
+    _load_tile(xq_ptr, xs_ptr, ids_ptr, eid_ptr, w_desc, ws_ptr, bufs, bars, pid_m, pid_n, N, nk, 0, 0, num_valid,
+               stride_xq, ks4, n_rows, x_div, K_TAIL, GATED)
+
+
+@gluon.jit
+def _mma_partition(bufs, bars, nk, mma: gl.constexpr, wsl: gl.constexpr, xsl: gl.constexpr, K_TAIL: gl.constexpr):
+    x_bufs, w_bufs, sw_bufs, sx_bufs = bufs
+    ready, empty, s_ready, s_empty = bars
+    S: gl.constexpr = x_bufs.type.shape[0]
+    BM: gl.constexpr = x_bufs.type.shape[1]
+    KB: gl.constexpr = x_bufs.type.shape[2]
+    BN: gl.constexpr = w_bufs.type.shape[1]
+    SB: gl.constexpr = KB // 8
+    SPC: gl.constexpr = 128 // KB
+    al: gl.constexpr = gl.DotOperandLayout(0, mma, 4)
+    bl: gl.constexpr = gl.DotOperandLayout(1, mma, 4)
+    acc = gl.zeros([BN, BM], gl.float32, mma)
+    # with K % 256 != 0 the last scale block covers fewer than SPC chunks
+    for jb in range(gl.cdiv(nk, SPC)):
+        ss = jb % 2
+        mbarrier.wait(s_ready.index(ss), jb // 2 & 1)
+        for c in gl.static_range(SPC):
+            kk = jb * SPC + c
+            if not K_TAIL or kk < nk:
+                st = kk % S
+                mbarrier.wait(ready.index(st), kk // S & 1)
+                a = w_bufs.index(st).load(al)
+                b = x_bufs.index(st).permute((1, 0)).load(bl)
+                sw = sw_bufs.index(ss).slice(c * SB, SB, dim=1).load(wsl).to(gl.float8e4nv, bitcast=True)
+                sx = sx_bufs.index(ss).slice(c * SB, SB, dim=1).load(xsl).to(gl.float8e4nv, bitcast=True)
+                mbarrier.arrive(empty.index(st))
+                acc = mma_nvf4(a, sw, b, sx, acc)
+        mbarrier.arrive(s_empty.index(ss))
+    return acc
+
+
+@gluon.jit
+def _e2m1x2(hi, lo):
+    # cvt packs its first source into the high nibble
+    return gl.inline_asm_elementwise("{ .reg .b8 t; cvt.rn.satfinite.e2m1x2.f32 t, $1, $2; cvt.u16.u8 $0, t; }",
+                                     "=h,r,r", [hi, lo], gl.uint16, True, 1)
+
+
+@gluon.jit
+def _gated_epilogue(acc, slot, pid_m, pid_n, wg_ptr, g_in, ids_ptr, q_ptr, qs_ptr, g_out, N, num_valid, stride_ge,
+                    ACT: gl.constexpr, act_alpha, act_limit):
+    """act(gate, up) of the tile's BN / 2 columns as NVFP4 under the global g_out: q [routes, N / 4], qs [routes, N / 32]."""
+    BN: gl.constexpr = acc.type.shape[0]
+    BM: gl.constexpr = acc.type.shape[1]
+    H2: gl.constexpr = BN // 2
+    G: gl.constexpr = H2 // 16
+    NW: gl.constexpr = gl.num_warps()
+    out_l: gl.constexpr = acc.type.layout
+    r = gl.arange(0, BN, layout=gl.SliceLayout(1, out_l))
+    col = gl.where(r < H2, pid_n * H2 + r, N // 2 + pid_n * H2 + r - H2)
+    wg = gl.load(wg_ptr + slot.to(gl.int64) * stride_ge + col).to(gl.float32) * g_in
+    # bf16 like the W4A16 gate_up output; n-contiguous 16-byte runs, swizzled so the lanes of a warp spread over banks
+    h_l: gl.constexpr = gl.SwizzledSharedLayout(8, 1, 8, [0, 1])
+    h_s = gl.allocate_shared_memory(gl.bfloat16, [BN, BM], h_l, (acc * wg[:, None]).to(gl.bfloat16))
+    # one thread holds 16 consecutive columns of one route
+    TG: gl.constexpr = min(32, G)
+    ql2: gl.constexpr = gl.BlockedLayout([16, 1], [TG, 32 // TG], [1, NW], [0, 1])
+    l4: gl.constexpr = gl.BlockedLayout([1, 8, 2, 1], [TG, 1, 1, 32 // TG], [1, 1, 1, NW], [2, 1, 0, 3])
+    gate = gl.convert_layout(gl.reshape(h_s.slice(0, H2).load(ql2), [G, 8, 2, BM]), l4, assert_trivial=True).to(gl.float32)
+    up = gl.convert_layout(gl.reshape(h_s.slice(H2, H2).load(ql2), [G, 8, 2, BM]), l4, assert_trivial=True).to(gl.float32)
+    if ACT == 1:
+        inner = 0.7978845608028654 * (gate + 0.044715 * gate * gate * gate)
+        # tanh(y) = 1 - 2 / (exp(2y) + 1)
+        h = 0.5 * gate * (2.0 - 2.0 / (gl.exp(2.0 * inner) + 1.0)) * up
+    elif ACT == 2 or ACT == 3:
+        gate = gl.minimum(gate, act_limit)
+        up = gl.minimum(gl.maximum(up, -act_limit), act_limit)
+        if ACT == 2:
+            up = up + 1.0
+        h = gate / (1.0 + gl.exp(-gate * act_alpha)) * up
+    else:
+        h = gate / (1.0 + gl.exp(-gate)) * up
+    h = h.to(gl.bfloat16).to(gl.float32)
+    amax = gl.max(gl.max(gl.abs(h), axis=2), axis=1)
+    sc = (amax / (6.0 * g_out)).to(gl.float8e4nv)
+    sf = sc.to(gl.float32) * g_out
+    inv = gl.where(sf > 0, 1.0 / sf, 0.0)
+    v = h * gl.expand_dims(gl.expand_dims(inv, 1), 2)
+    lo, hi = gl.split(gl.permute(v, (0, 1, 3, 2)))
+    q = _e2m1x2(hi, lo).to(gl.uint8)
+    ql: gl.constexpr = q.type.layout
+    g_i = gl.arange(0, G, layout=gl.SliceLayout(1, gl.SliceLayout(2, ql)))
+    b_i = gl.arange(0, 8, layout=gl.SliceLayout(0, gl.SliceLayout(2, ql)))
+    rq = gl.load(ids_ptr + pid_m * BM + gl.arange(0, BM, layout=gl.SliceLayout(0, gl.SliceLayout(1, ql))))
+    q_cols = pid_n * (H2 // 2) + g_i[:, None] * 8 + b_i[None, :]
+    q_ptrs = q_ptr + (rq.to(gl.int64) * (N // 4))[None, None, :] + q_cols[:, :, None]
+    gl.store(q_ptrs, q, mask=(rq < num_valid)[None, None, :])
+    sl: gl.constexpr = sc.type.layout
+    sg = gl.arange(0, G, layout=gl.SliceLayout(1, sl))
+    rs = gl.load(ids_ptr + pid_m * BM + gl.arange(0, BM, layout=gl.SliceLayout(0, sl)))
+    s_ptrs = qs_ptr + (rs.to(gl.int64) * (N // 32))[None, :] + (pid_n * G + sg)[:, None]
+    gl.store(s_ptrs, sc.to(gl.uint8, bitcast=True), mask=(rs < num_valid)[None, :])
+
+
+@gluon.jit
+def _w4a4_moe_kernel(
+    xq_ptr, xs_ptr, g_in_ptr, w_desc, ws_ptr, wg_ptr, out_ptr, topk_weights_ptr,
+    sorted_token_ids_ptr, expert_ids_ptr, num_tokens_post_padded_ptr,
+    q_ptr, qs_ptr, g_out_ptr, act_alpha, act_limit,
+    N, K, EM, num_valid_tokens, ks4, n_rows,
+    stride_xq, stride_ge, stride_out, stride_tw,
+    BLOCK_SIZE_M: gl.constexpr, BLOCK_SIZE_N: gl.constexpr, BLOCK_SIZE_K: gl.constexpr,
+    GROUP_SIZE_M: gl.constexpr, N_GROUP: gl.constexpr, STAGES: gl.constexpr, MMA_WARPS_N: gl.constexpr,
+    LOAD_WARPS: gl.constexpr, LOAD_REGS: gl.constexpr, K_TAIL: gl.constexpr, GATED: gl.constexpr,
+    ACT: gl.constexpr, top_k: gl.constexpr, wsl: gl.constexpr, xsl: gl.constexpr,
+):
+    BM: gl.constexpr = BLOCK_SIZE_M
+    BN: gl.constexpr = BLOCK_SIZE_N
+    KB: gl.constexpr = BLOCK_SIZE_K // 2
+    num_pid_m = gl.cdiv(EM, BM)
+    pid_m, pid_n = _tile(gl.program_id(0), num_pid_m, gl.cdiv(N, BN), GROUP_SIZE_M, N_GROUP)
+    if pid_m >= gl.cdiv(gl.load(num_tokens_post_padded_ptr), BM):
+        return
+    nk = K // BLOCK_SIZE_K
+    # gate_up reads one x row per token, down one row per route
+    x_div = top_k if GATED else 1
+
+    x_bufs = gl.allocate_shared_memory(gl.uint8, [STAGES, BM, KB], gl.NVMMASharedLayout.get_default_for([BM, KB], gl.uint8))
+    w_bufs = gl.allocate_shared_memory(gl.uint8, [STAGES, BN, KB], w_desc.layout)
+    sw_bufs = gl.allocate_shared_memory(gl.uint8, [2, BN, 16], SCALE_SMEM_LAYOUT)
+    sx_bufs = gl.allocate_shared_memory(gl.uint8, [2, BM, 16], SCALE_SMEM_LAYOUT)
+    ready = gl.allocate_shared_memory(gl.int64, [STAGES, 1], mbarrier.MBarrierLayout())
+    empty = gl.allocate_shared_memory(gl.int64, [STAGES, 1], mbarrier.MBarrierLayout())
+    s_ready = gl.allocate_shared_memory(gl.int64, [2, 1], mbarrier.MBarrierLayout())
+    s_empty = gl.allocate_shared_memory(gl.int64, [2, 1], mbarrier.MBarrierLayout())
+    for i in gl.static_range(STAGES):
+        mbarrier.init(ready.index(i), count=1)
+        mbarrier.init(empty.index(i), count=1)
+    for i in gl.static_range(2):
+        mbarrier.init(s_ready.index(i), count=1)
+        mbarrier.init(s_empty.index(i), count=1)
+
+    NW: gl.constexpr = gl.num_warps()
+    mma: gl.constexpr = gl.NVMMADistributedLayout(version=[2, 0], warps_per_cta=[MMA_WARPS_N, NW // MMA_WARPS_N],
+                                                  instr_shape=[16, 8])
+    # tuples go inline: a tuple stored in a variable turns its constexprs into runtime values
+    acc = gl.warp_specialize([
+        (_mma_partition, ((x_bufs, w_bufs, sw_bufs, sx_bufs), (ready, empty, s_ready, s_empty), nk, mma, wsl, xsl, K_TAIL)),
+        (_load_partition, (xq_ptr, xs_ptr, sorted_token_ids_ptr, expert_ids_ptr, w_desc, ws_ptr,
+                           (x_bufs, w_bufs, sw_bufs, sx_bufs), (ready, empty, s_ready, s_empty), N, nk, num_pid_m,
+                           num_valid_tokens, stride_xq, ks4, n_rows, x_div, K_TAIL, GROUP_SIZE_M, GATED, N_GROUP)),
+    ], [LOAD_WARPS], [LOAD_REGS])
+
+    out_l: gl.constexpr = gl.BlockedLayout([8, 1], [BN // 8, 32 * 8 // BN], [1, NW], [0, 1])
+    acc = gl.convert_layout(acc, out_l)
+    slot = gl.load(expert_ids_ptr + pid_m)
+    g_in = gl.load(g_in_ptr)
+    if GATED:
+        _gated_epilogue(acc, slot, pid_m, pid_n, wg_ptr, g_in, sorted_token_ids_ptr, q_ptr, qs_ptr, gl.load(g_out_ptr),
+                        N, num_valid_tokens, stride_ge, ACT, act_alpha, act_limit)
+    else:
+        offs_n = pid_n * BN + gl.arange(0, BN, layout=gl.SliceLayout(1, out_l))
+        n_mask = offs_n < N
+        route = gl.load(sorted_token_ids_ptr + pid_m * BM + gl.arange(0, BM, layout=gl.SliceLayout(0, out_l)))
+        route_mask = route < num_valid_tokens
+        wg = gl.load(wg_ptr + slot.to(gl.int64) * stride_ge + offs_n, mask=n_mask, other=0.0).to(gl.float32) * g_in
+        tw = gl.load(topk_weights_ptr + route * stride_tw, mask=route_mask, other=0.0)
+        acc = acc * wg[:, None] * tw[None, :]
+        out_ptrs = out_ptr + (route // top_k).to(gl.int64)[None, :] * stride_out + offs_n[:, None]
+        gl.atomic_add(out_ptrs, acc.to(out_ptr.dtype.element_ty), mask=route_mask[None, :] & n_mask[:, None], sem="relaxed")
+
+
+# Persistent down GEMM: the 256 x 128 tile and 8 MMA warps of the kernel above, with the routing-weight scale and the atomics on
+# separate epilogue warps. The accumulator is two 128 x 128 halves that go to the epilogue one at a time through one 32 KB buffer,
+# so the next tile's mainloop overlaps the atomics of the last half.
+
+
+@gluon.jit
+def _num_tiles(ntpp_ptr, N, BM: gl.constexpr, BN: gl.constexpr):
+    num_pid_m = gl.cdiv(gl.load(ntpp_ptr), BM)
+    num_pid_n = gl.cdiv(N, BN)
+    return num_pid_m, num_pid_n, num_pid_m * num_pid_n
+
+
+@gluon.jit
+def _down_load_partition(xq_ptr, xs_ptr, ids_ptr, eid_ptr, ntpp_ptr, w_desc, ws_ptr, bufs, bars, N, nk, num_valid,
+                         stride_xq, ks4, n_rows, K_TAIL: gl.constexpr, N_GROUP: gl.constexpr, GROUP_SIZE_M: gl.constexpr):
+    BM: gl.constexpr = bufs[0].type.shape[1]
+    KB: gl.constexpr = bufs[0].type.shape[2]
+    BN: gl.constexpr = bufs[1].type.shape[1]
+    SPC: gl.constexpr = 128 // KB
+    num_pid_m, num_pid_n, num_tiles = _num_tiles(ntpp_ptr, N, BM, BN)
+    # it / jt count chunks and scale blocks across tiles, so the rings keep their phase from one tile to the next
+    it = 0
+    jt = 0
+    for t in range(gl.program_id(0), num_tiles, gl.num_programs(0)):
+        pid_m, pid_n = _tile(t, num_pid_m, num_pid_n, GROUP_SIZE_M, N_GROUP)
+        _load_tile(xq_ptr, xs_ptr, ids_ptr, eid_ptr, w_desc, ws_ptr, (bufs[0], bufs[1], bufs[2], bufs[3]),
+                   (bars[0], bars[1], bars[2], bars[3]), pid_m, pid_n, N, nk, it, jt, num_valid, stride_xq, ks4,
+                   n_rows, 1, K_TAIL, False)
+        it += nk
+        jt += gl.cdiv(nk, SPC)
+
+
+@gluon.jit
+def _down_epi_partition(ids_ptr, ntpp_ptr, out_ptr, bufs, bars, N, num_valid, stride_out, top_k: gl.constexpr,
+                        N_GROUP: gl.constexpr, GROUP_SIZE_M: gl.constexpr):
+    h_buf = bufs[4]
+    acc_ready, acc_empty = bars[4], bars[5]
+    HB: gl.constexpr = h_buf.type.shape[0]
+    BM: gl.constexpr = h_buf.type.shape[1]
+    BN: gl.constexpr = 2 * HB
+    CH: gl.constexpr = 64
+    NWE: gl.constexpr = gl.num_warps()
+    # 8 consecutive n of one route per thread: one 16-byte bf16 atomic; 64 rows per pass keep the registers low
+    el: gl.constexpr = gl.BlockedLayout([8, 1], [CH // 8, 32 * 8 // CH], [1, NWE], [0, 1])
+    offs_cl = gl.arange(0, CH, layout=gl.SliceLayout(1, el))
+    offs_tl = gl.arange(0, BM, layout=gl.SliceLayout(0, el))
+    num_pid_m, num_pid_n, num_tiles = _num_tiles(ntpp_ptr, N, BM, BN)
+    hc = 0
+    for t in range(gl.program_id(0), num_tiles, gl.num_programs(0)):
+        pid_m, pid_n = _tile(t, num_pid_m, num_pid_n, GROUP_SIZE_M, N_GROUP)
+        route = gl.load(ids_ptr + pid_m * BM + offs_tl)
+        rmask = (route < num_valid)[None, :]
+        rows = (route // top_k).to(gl.int64)[None, :] * stride_out
+        for half in gl.static_range(2):
+            mbarrier.wait(acc_ready.index(0), hc & 1)
+            for ch in gl.static_range(HB // CH):
+                offs_n = pid_n * BN + half * HB + ch * CH + offs_cl
+                h = h_buf.slice(ch * CH, CH).load(el)
+                gl.atomic_add(out_ptr + rows + offs_n[:, None], h, mask=rmask & (offs_n < N)[:, None], sem="relaxed")
+            mbarrier.arrive(acc_empty.index(0))
+            hc += 1
+
+
+@gluon.jit
+def _down_mma_partition(bufs, bars, ids_ptr, eid_ptr, ntpp_ptr, wg_ptr, g_in_ptr, tw_ptr, N, nk, num_valid, stride_ge,
+                        stride_tw, mma: gl.constexpr, wsl: gl.constexpr, xsl: gl.constexpr, K_TAIL: gl.constexpr,
+                        N_GROUP: gl.constexpr, GROUP_SIZE_M: gl.constexpr):
+    x_bufs, w_bufs, sw_bufs, sx_bufs, h_buf = bufs
+    ready, empty, s_ready, s_empty, acc_ready, acc_empty = bars
+    S: gl.constexpr = x_bufs.type.shape[0]
+    BM: gl.constexpr = x_bufs.type.shape[1]
+    KB: gl.constexpr = x_bufs.type.shape[2]
+    BN: gl.constexpr = w_bufs.type.shape[1]
+    HB: gl.constexpr = BN // 2
+    SB: gl.constexpr = KB // 8
+    SPC: gl.constexpr = 128 // KB
+    al: gl.constexpr = gl.DotOperandLayout(0, mma, 4)
+    bl: gl.constexpr = gl.DotOperandLayout(1, mma, 4)
+    offs_nl = gl.arange(0, HB, layout=gl.SliceLayout(1, mma))
+    offs_tl = gl.arange(0, BM, layout=gl.SliceLayout(0, mma))
+    g_in = gl.load(g_in_ptr)
+    num_pid_m, num_pid_n, num_tiles = _num_tiles(ntpp_ptr, N, BM, BN)
+    it = 0
+    jt = 0
+    hc = 0
+    for t in range(gl.program_id(0), num_tiles, gl.num_programs(0)):
+        pid_m, pid_n = _tile(t, num_pid_m, num_pid_n, GROUP_SIZE_M, N_GROUP)
+        acc0 = gl.zeros([HB, BM], gl.float32, mma)
+        acc1 = gl.zeros([HB, BM], gl.float32, mma)
+        for jb in range(gl.cdiv(nk, SPC)):
+            j = jt + jb
+            ss = j % 2
+            mbarrier.wait(s_ready.index(ss), j // 2 & 1)
+            for c in gl.static_range(SPC):
+                kk = jb * SPC + c
+                if not K_TAIL or kk < nk:
+                    i = it + kk
+                    st = i % S
+                    mbarrier.wait(ready.index(st), i // S & 1)
+                    b = x_bufs.index(st).permute((1, 0)).load(bl)
+                    sx = sx_bufs.index(ss).slice(c * SB, SB, dim=1).load(xsl).to(gl.float8e4nv, bitcast=True)
+                    a0 = w_bufs.index(st).slice(0, HB).load(al)
+                    a1 = w_bufs.index(st).slice(HB, HB).load(al)
+                    sws = sw_bufs.index(ss).slice(c * SB, SB, dim=1)
+                    sw0 = sws.slice(0, HB).load(wsl).to(gl.float8e4nv, bitcast=True)
+                    sw1 = sws.slice(HB, HB).load(wsl).to(gl.float8e4nv, bitcast=True)
+                    mbarrier.arrive(empty.index(st))
+                    acc0 = mma_nvf4(a0, sw0, b, sx, acc0)
+                    acc1 = mma_nvf4(a1, sw1, b, sx, acc1)
+            mbarrier.arrive(s_empty.index(ss))
+        it += nk
+        jt += gl.cdiv(nk, SPC)
+        slot = gl.load(eid_ptr + pid_m)
+        route = gl.load(ids_ptr + pid_m * BM + offs_tl)
+        rmask = route < num_valid
+        xs_ = gl.load(tw_ptr + route * stride_tw, mask=rmask, other=0.0) * g_in
+        wrow = wg_ptr + slot.to(gl.int64) * stride_ge
+        n0 = pid_n * BN + offs_nl
+        wg0 = gl.load(wrow + n0, mask=n0 < N, other=0.0).to(gl.float32)
+        h0 = (acc0 * wg0[:, None] * xs_[None, :]).to(gl.bfloat16)
+        mbarrier.wait(acc_empty.index(0), (hc & 1) ^ 1)
+        h_buf.store(h0)
+        mbarrier.arrive(acc_ready.index(0))
+        hc += 1
+        n1 = n0 + HB
+        wg1 = gl.load(wrow + n1, mask=n1 < N, other=0.0).to(gl.float32)
+        h1 = (acc1 * wg1[:, None] * xs_[None, :]).to(gl.bfloat16)
+        mbarrier.wait(acc_empty.index(0), (hc & 1) ^ 1)
+        h_buf.store(h1)
+        mbarrier.arrive(acc_ready.index(0))
+        hc += 1
+
+
+@gluon.jit
+def _w4a4_down_persistent_kernel(
+    xq_ptr, xs_ptr, g_in_ptr, w_desc, ws_ptr, wg_ptr, out_ptr, topk_weights_ptr,
+    sorted_token_ids_ptr, expert_ids_ptr, num_tokens_post_padded_ptr,
+    N, K, num_valid_tokens, ks4, n_rows, stride_xq, stride_ge, stride_out, stride_tw,
+    BLOCK_SIZE_M: gl.constexpr, BLOCK_SIZE_N: gl.constexpr, BLOCK_SIZE_K: gl.constexpr, STAGES: gl.constexpr,
+    K_TAIL: gl.constexpr, N_GROUP: gl.constexpr, GROUP_SIZE_M: gl.constexpr, MMA_WARPS_N: gl.constexpr,
+    LOAD_WARPS: gl.constexpr, LOAD_REGS: gl.constexpr, EPI_WARPS: gl.constexpr, EPI_REGS: gl.constexpr,
+    top_k: gl.constexpr, wsl: gl.constexpr, xsl: gl.constexpr,
+):
+    BM: gl.constexpr = BLOCK_SIZE_M
+    BN: gl.constexpr = BLOCK_SIZE_N
+    KB: gl.constexpr = BLOCK_SIZE_K // 2
+    nk = K // BLOCK_SIZE_K
+    x_bufs = gl.allocate_shared_memory(gl.uint8, [STAGES, BM, KB], gl.NVMMASharedLayout.get_default_for([BM, KB], gl.uint8))
+    w_bufs = gl.allocate_shared_memory(gl.uint8, [STAGES, BN, KB], w_desc.layout)
+    sw_bufs = gl.allocate_shared_memory(gl.uint8, [2, BN, 16], SCALE_SMEM_LAYOUT)
+    sx_bufs = gl.allocate_shared_memory(gl.uint8, [2, BM, 16], SCALE_SMEM_LAYOUT)
+    # n-contiguous and swizzled, so the epilogue's 8-row reads and the MMA-layout stores spread over banks
+    h_buf = gl.allocate_shared_memory(gl.bfloat16, [BN // 2, BM], gl.SwizzledSharedLayout(8, 1, 8, [0, 1]))
+    ready = gl.allocate_shared_memory(gl.int64, [STAGES, 1], mbarrier.MBarrierLayout())
+    empty = gl.allocate_shared_memory(gl.int64, [STAGES, 1], mbarrier.MBarrierLayout())
+    s_ready = gl.allocate_shared_memory(gl.int64, [2, 1], mbarrier.MBarrierLayout())
+    s_empty = gl.allocate_shared_memory(gl.int64, [2, 1], mbarrier.MBarrierLayout())
+    acc_ready = gl.allocate_shared_memory(gl.int64, [1, 1], mbarrier.MBarrierLayout())
+    acc_empty = gl.allocate_shared_memory(gl.int64, [1, 1], mbarrier.MBarrierLayout())
+    for i in gl.static_range(STAGES):
+        mbarrier.init(ready.index(i), count=1)
+        mbarrier.init(empty.index(i), count=1)
+    for i in gl.static_range(2):
+        mbarrier.init(s_ready.index(i), count=1)
+        mbarrier.init(s_empty.index(i), count=1)
+    mbarrier.init(acc_ready.index(0), count=1)
+    mbarrier.init(acc_empty.index(0), count=1)
+
+    NW: gl.constexpr = gl.num_warps()
+    mma: gl.constexpr = gl.NVMMADistributedLayout(version=[2, 0], warps_per_cta=[MMA_WARPS_N, NW // MMA_WARPS_N],
+                                                  instr_shape=[16, 8])
+    gl.warp_specialize([
+        (_down_mma_partition, ((x_bufs, w_bufs, sw_bufs, sx_bufs, h_buf), (ready, empty, s_ready, s_empty, acc_ready, acc_empty),
+                               sorted_token_ids_ptr, expert_ids_ptr, num_tokens_post_padded_ptr, wg_ptr, g_in_ptr,
+                               topk_weights_ptr, N, nk, num_valid_tokens, stride_ge, stride_tw, mma, wsl, xsl, K_TAIL,
+                               N_GROUP, GROUP_SIZE_M)),
+        (_down_load_partition, (xq_ptr, xs_ptr, sorted_token_ids_ptr, expert_ids_ptr, num_tokens_post_padded_ptr, w_desc,
+                                ws_ptr, (x_bufs, w_bufs, sw_bufs, sx_bufs, h_buf),
+                                (ready, empty, s_ready, s_empty, acc_ready, acc_empty), N, nk, num_valid_tokens,
+                                stride_xq, ks4, n_rows, K_TAIL, N_GROUP, GROUP_SIZE_M)),
+        (_down_epi_partition, (sorted_token_ids_ptr, num_tokens_post_padded_ptr, out_ptr,
+                               (x_bufs, w_bufs, sw_bufs, sx_bufs, h_buf), (ready, empty, s_ready, s_empty, acc_ready, acc_empty),
+                               N, num_valid_tokens, stride_out, top_k, N_GROUP, GROUP_SIZE_M)),
+    ], [LOAD_WARPS, EPI_WARPS], [LOAD_REGS, EPI_REGS])
+
+
+def supports(H: int, I: int) -> bool:
+    """Whether the GEMMs take a layer of this hidden / intermediate size."""
+    # K is a whole number of 64-wide k blocks for both GEMMs; the gated tile pairs 64 gate rows with 64 up rows at least
+    return H % 64 == 0 and I % 64 == 0
+
+
+def _weight_desc(packed, box_n, BK):
+    S, N, KB = packed.shape
+    layout = gl.NVMMASharedLayout.get_default_for([box_n, BK // 2], gl.uint8)
+    return TensorDescriptor.from_tensor(packed.view(S * N, KB), [box_n, BK // 2], layout)
+
+
+def _scale_layouts(cfg, w_rows):
+    nw, wn = cfg["num_warps"], cfg["MMA_WARPS_N"]
+    kg = cfg["BLOCK_SIZE_K"] // 16
+    return scale_layout(0, w_rows, kg, [wn, nw // wn]), scale_layout(1, cfg["BLOCK_SIZE_M"], kg, [wn, nw // wn])
+
+
+def gate_up_gemm(xq, xs, g_in, packed, scale, glob, q, qs, g_out, sorted_ids, expert_ids, num_tokens_post_padded,
+                 num_valid_tokens, top_k, activation, act_alpha, act_limit, cfg):
+    """q, qs = NVFP4 under g_out of act(gate, up), where [gate; up] = (xq * xs * g_in)[route // top_k] W4A4 packed[expert]^T * glob."""
+    S, N, KB = packed.shape
+    K = KB * 2
+    BM, BN, BK = cfg["BLOCK_SIZE_M"], cfg["BLOCK_SIZE_N"], cfg["BLOCK_SIZE_K"]
+    assert K % BK == 0 and (N // 2) % (BN // 2) == 0 and xs.stride(0) == K // 16
+    assert q.stride(0) == N // 4 and qs.stride(0) == N // 32
+    wsl, xsl = _scale_layouts(cfg, BN)
+    grid = triton.cdiv(sorted_ids.shape[0], BM) * triton.cdiv(N, BN)
+    _w4a4_moe_kernel[(grid,)](
+        xq, xs.view(torch.uint8), g_in, _weight_desc(packed, BN // 2, BK), scale.view(torch.uint8), glob, None, None,
+        sorted_ids, expert_ids, num_tokens_post_padded, q, qs.view(torch.uint8), g_out, act_alpha, act_limit,
+        N, K, sorted_ids.shape[0], num_valid_tokens, K // 64, S * N, xq.stride(0), glob.stride(0), 0, 0,
+        K_TAIL=K % 256 != 0, GATED=True, ACT=ACTIVATIONS[activation], top_k=top_k, wsl=wsl, xsl=xsl, **cfg,
+    )
+
+
+def down_gemm(xq, xs, g_in, packed, scale, glob, out, topk_weights_flat, sorted_ids, expert_ids, num_tokens_post_padded,
+              num_valid_tokens, top_k, cfg):
+    """out[route // top_k] += (xq * xs * g_in)[route] W4A4 packed[expert]^T * glob * topk_weights_flat[route]; out starts zeroed."""
+    S, N, KB = packed.shape
+    K = KB * 2
+    cfg = dict(cfg)
+    persistent = cfg.pop("PERSISTENT")
+    BM, BN, BK = cfg["BLOCK_SIZE_M"], cfg["BLOCK_SIZE_N"], cfg["BLOCK_SIZE_K"]
+    assert K % BK == 0 and xs.stride(0) == K // 16
+    args = (xq, xs.view(torch.uint8), g_in, _weight_desc(packed, BN, BK), scale.view(torch.uint8), glob, out,
+            topk_weights_flat, sorted_ids, expert_ids, num_tokens_post_padded)
+    if persistent:
+        wsl, xsl = _scale_layouts(cfg, BN // 2)
+        grid = torch.cuda.get_device_properties(out.device).multi_processor_count
+        _w4a4_down_persistent_kernel[(grid,)](
+            *args, N, K, num_valid_tokens, K // 64, S * N, xq.stride(0), glob.stride(0), out.stride(0),
+            topk_weights_flat.stride(0), K_TAIL=K % 256 != 0, top_k=top_k, wsl=wsl, xsl=xsl, **cfg,
+        )
+        return
+    wsl, xsl = _scale_layouts(cfg, BN)
+    grid = triton.cdiv(sorted_ids.shape[0], BM) * triton.cdiv(N, BN)
+    _w4a4_moe_kernel[(grid,)](
+        *args, None, None, None, 1.0, 1.0,
+        N, K, sorted_ids.shape[0], num_valid_tokens, K // 64, S * N, xq.stride(0), glob.stride(0), out.stride(0),
+        topk_weights_flat.stride(0), K_TAIL=K % 256 != 0, GATED=False, ACT=0, top_k=top_k, wsl=wsl, xsl=xsl, **cfg,
+    )
diff --git a/python/freetoken/kernel/gluon/nvfp4/utils.py b/python/freetoken/kernel/gluon/nvfp4/utils.py
index 430c146d4..11426ec56 100644
--- a/python/freetoken/kernel/gluon/nvfp4/utils.py
+++ b/python/freetoken/kernel/gluon/nvfp4/utils.py
@@ -2,6 +2,11 @@
 
 from triton.experimental import gluon
 from triton.experimental.gluon import language as gl
+from triton.experimental.gluon.language import _core as ttgl
+from triton.experimental.gluon.language._layouts import (
+    DotOperandLayout,
+    NVMMADistributedLayout,
+)
 
 
 @gluon.jit
@@ -41,3 +46,37 @@ def e2m1x8_bf16x2(q, s2):
     w2 = bf16x2_mul(prmt(o3, o2, 0x7632), s2)
     w3 = bf16x2_mul(prmt(o1, o0, 0x7632), s2)
     return w0, w1, w2, w3
+
+
+def scale_layout(op_idx: int, rows: int, kgroups: int, mma_warps: list[int]) -> gl.DistributedLinearLayout:
+    """Layout of a [rows, K / 16] e4m3 scale tile for operand op_idx of an NVFP4 mma_v2; mirrors Triton's getSM120DotScaledScaleLayout."""
+    def bits(n):
+        return n.bit_length() - 1
+
+    w0, w1 = mma_warps
+    reg = [[0, 1 << i] for i in range(bits(kgroups))]
+    if op_idx == 0:
+        lane = [[8, 0], [0, 0], [1, 0], [2, 0], [4, 0]]
+        warp = [[0, 0]] * bits(w1) + [[16 << i, 0] for i in range(bits(w0))]
+        cover = 16 * w0
+    else:
+        lane = [[0, 0], [0, 0], [1, 0], [2, 0], [4, 0]]
+        warp = [[8 << i, 0] for i in range(bits(w1))] + [[0, 0]] * bits(w0)
+        cover = 8 * w1
+    reg += [[cover << i, 0] for i in range(max(0, bits(rows) - bits(cover)))]
+    # bases past the tile's rows broadcast
+    lane = [[b[0] if b[0] < rows else 0, b[1]] for b in lane]
+    warp = [[b[0] if b[0] < rows else 0, b[1]] for b in warp]
+    return gl.DistributedLinearLayout(reg, lane, warp, [], [rows, kgroups])
+
+
+@ttgl.builtin
+def mma_nvf4(a, a_scale, b, b_scale, acc, _semantic=None):
+    """acc += (a * a_scale) @ (b * b_scale) on sm_120 (mma kind::mxf4nvf4, scale_vec::4X).
+
+    a [M, K/2] and b [K/2, N] are packed e2m1 in DotOperandLayout(k_width=4) of acc's mma_v2 layout; the e4m3 scales are [M, K/16] / [N, K/16] in scale_layout.
+    """
+    assert isinstance(acc.type.layout, NVMMADistributedLayout) and acc.type.layout.version == [2, 0]
+    assert isinstance(a.type.layout, DotOperandLayout) and isinstance(b.type.layout, DotOperandLayout)
+    out = _semantic.dot_scaled(a, a_scale, "e2m1", b, b_scale, "e2m1", acc, False, True, True, ttgl.float32)
+    return ttgl.tensor(out.handle, acc.type)
diff --git a/python/freetoken/kernel/nvfp4_moe.py b/python/freetoken/kernel/nvfp4_moe.py
index 0aa494222..df8941072 100644
--- a/python/freetoken/kernel/nvfp4_moe.py
+++ b/python/freetoken/kernel/nvfp4_moe.py
@@ -1,6 +1,6 @@
-"""Kernel dispatch for the inline-dequant NVFP4 MoE GEMMs: which kernel runs, with which tiles and switches, and its launch.
+"""Kernel dispatch for the NVFP4 MoE GEMMs: which kernel runs, with which tiles and switches, and its launch.
 
-The kernels are templates (:mod:`freetoken.kernel.triton.nvfp4_fused_moe`, :mod:`freetoken.kernel.gluon.nvfp4.moe.sm120`);
+The kernels are templates (:mod:`freetoken.kernel.triton.nvfp4_fused_moe`, :mod:`freetoken.kernel.gluon.nvfp4.moe.sm120`, :mod:`freetoken.kernel.gluon.nvfp4.moe.sm120_w4a4`);
 :mod:`freetoken.moe.fused_nvfp4` strings the GEMMs into the MoE layer.
 """
 
@@ -388,3 +388,73 @@ def prefill_gemm(
         SCALE_WORDS=scale_words,
         **cfg,
     )
+
+
+@functools.lru_cache(maxsize=None)
+def _w4a4_prefill(device_index: int):
+    """The Gluon module of the W4A4 prefill GEMMs on this GPU, or None."""
+    from freetoken.kernel import gluon
+
+    if torch.cuda.get_device_capability(device_index) != (12, 0) or not gluon.available():
+        return None
+    from freetoken.kernel.gluon.nvfp4.moe import sm120_w4a4
+
+    return sm120_w4a4
+
+
+def w4a4_available(device_index: int) -> bool:
+    return _w4a4_prefill(device_index) is not None
+
+
+# Below this many tokens per expert both prefills stream the weights at the same speed (RTX PRO 6000 and RTX 5090 sweeps over the
+# NVFP4 MoE shapes); W4A16 keeps the exact activations there
+_W4A4_MIN_TOKENS_PER_EXPERT = 4
+# (min rho, tokens per expert where BLOCK_SIZE_M becomes 128), from the same sweeps
+_W4A4_BLOCK_M_STEPS = ((3.0, 128), (0.0, 64))
+
+
+@functools.lru_cache(maxsize=None)
+def _w4a4_device(device_index: int) -> Tuple[float, int]:
+    props = torch.cuda.get_device_properties(device_index)
+    rho = _tensor_rho(props.name, props.major, props.minor)
+    return next(t for min_rho, t in _W4A4_BLOCK_M_STEPS if rho >= min_rho), props.L2_cache_size
+
+
+def w4a4_prefill_config(M: int, top_k: int, num_experts: int, H: int, I: int, activation: str, dtype: torch.dtype,
+                        device_index: int) -> Dict[str, Any] | None:
+    """Tiles of the W4A4 prefill GEMMs, or None where the W4A16 prefill runs instead."""
+    w4a4 = _w4a4_prefill(device_index)
+    if w4a4 is None or dtype != torch.bfloat16 or activation not in w4a4.ACTIVATIONS or not w4a4.supports(H, I):
+        return None
+    to128, l2_bytes = _w4a4_device(device_index)
+    per_expert = M * top_k / num_experts
+    if per_expert < _W4A4_MIN_TOKENS_PER_EXPERT:
+        return None
+    block_m = 128 if per_expert >= to128 else 64
+    common = dict(BLOCK_SIZE_M=block_m, GROUP_SIZE_M=8, MMA_WARPS_N=4, LOAD_WARPS=4, num_warps=8)
+    # the gated tile holds BN / 2 gate rows and their up rows, so I must split into whole half tiles
+    gate_up = dict(common, BLOCK_SIZE_N=256 if I % 128 == 0 else 128, BLOCK_SIZE_K=128 if H % 128 == 0 else 64,
+                   N_GROUP=0, STAGES=3, LOAD_REGS=88)
+    if I < 1024:
+        # short K: the persistent kernel overlaps a tile's atomics with the next mainloop; N strips keep the [M, H] output
+        # in L2 once it outgrows half of it
+        down = dict(common, BLOCK_SIZE_N=256, BLOCK_SIZE_K=64, STAGES=4, N_GROUP=2 if M * H * 2 > l2_bytes // 2 else 0,
+                    LOAD_REGS=56, EPI_WARPS=4, EPI_REGS=72, PERSISTENT=True)
+    else:
+        down = dict(common, BLOCK_SIZE_N=256 if H % 256 == 0 else 128, BLOCK_SIZE_K=128 if I % 128 == 0 else 64,
+                    N_GROUP=0, STAGES=3, LOAD_REGS=88, PERSISTENT=False)
+    return dict(BLOCK_SIZE_M=block_m, gate_up=gate_up, down=down)
+
+
+def w4a4_gate_up_gemm(xq, xs, g_in, packed, scale, glob, q, qs, g_out, sorted_ids, expert_ids, num_tokens_post_padded,
+                      num_valid_tokens, top_k, activation, act_alpha, act_limit, cfg):
+    _w4a4_prefill(xq.device.index).gate_up_gemm(
+        xq, xs, g_in, packed, scale, glob, q, qs, g_out, sorted_ids, expert_ids, num_tokens_post_padded,
+        num_valid_tokens, top_k, activation, act_alpha, act_limit, cfg)
+
+
+def w4a4_down_gemm(xq, xs, g_in, packed, scale, glob, out, topk_weights_flat, sorted_ids, expert_ids,
+                   num_tokens_post_padded, num_valid_tokens, top_k, cfg):
+    _w4a4_prefill(xq.device.index).down_gemm(
+        xq, xs, g_in, packed, scale, glob, out, topk_weights_flat, sorted_ids, expert_ids, num_tokens_post_padded,
+        num_valid_tokens, top_k, cfg)
diff --git a/python/freetoken/kernel/triton/nvfp4_act_quant.py b/python/freetoken/kernel/triton/nvfp4_act_quant.py
new file mode 100644
index 000000000..b52cfbe89
--- /dev/null
+++ b/python/freetoken/kernel/triton/nvfp4_act_quant.py
@@ -0,0 +1,42 @@
+"""NVFP4 activation quantization under a static global scale: x ~= e2m1 * e4m3 block scale (per 16) * g."""
+
+from __future__ import annotations
+
+import torch
+import triton
+import triton.language as tl
+
+
+@triton.jit
+def _e2m1x2(hi, lo):
+    # sm_100+ / sm_120a; cvt packs its first source into the high nibble
+    return tl.inline_asm_elementwise(
+        "{ .reg .b8 t; cvt.rn.satfinite.e2m1x2.f32 t, $1, $2; cvt.u16.u8 $0, t; }",
+        "=h,r,r", [hi, lo], dtype=tl.uint16, is_pure=True, pack=1)
+
+
+@triton.jit
+def _quant_kernel(x_ptr, q_ptr, s_ptr, g_ptr, K, stride_x, BLOCK_K: tl.constexpr):
+    row = tl.program_id(0).to(tl.int64)
+    offs = tl.arange(0, BLOCK_K)
+    x = tl.load(x_ptr + row * stride_x + offs, mask=offs < K, other=0.0).to(tl.float32)
+    g = tl.load(g_ptr)
+    xb = tl.reshape(x, [BLOCK_K // 16, 16])
+    s = (tl.max(tl.abs(xb), 1) / (6.0 * g)).to(tl.float8e4nv)
+    sf = s.to(tl.float32) * g
+    inv = tl.where(sf > 0, 1.0 / sf, 0.0)
+    lo, hi = tl.split(tl.reshape(xb * inv[:, None], [BLOCK_K // 2, 2]))
+    oq = tl.arange(0, BLOCK_K // 2)
+    tl.store(q_ptr + row * (K // 2) + oq, _e2m1x2(hi, lo).to(tl.uint8), mask=oq < K // 2)
+    os_ = tl.arange(0, BLOCK_K // 16)
+    tl.store(s_ptr + row * (K // 16) + os_, s.to(tl.uint8, bitcast=True), mask=os_ < K // 16)
+
+
+def quantize_nvfp4(x: torch.Tensor, g: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor]:
+    """x [R, K] -> (codes [R, K / 2] uint8, low nibble first; block scales [R, K / 16] e4m3) under the fp32 global g [1]; K % 16 == 0."""
+    R, K = x.shape
+    q = torch.empty(R, K // 2, dtype=torch.uint8, device=x.device)
+    s = torch.empty(R, K // 16, dtype=torch.uint8, device=x.device)
+    if R:
+        _quant_kernel[(R,)](x, q, s, g, K, x.stride(0), BLOCK_K=triton.next_power_of_2(K), num_warps=4)
+    return q, s.view(torch.float8_e4m3fn)
diff --git a/python/freetoken/layers/quantization/moe/nvfp4.py b/python/freetoken/layers/quantization/moe/nvfp4.py
index 1cd1a1ac7..0dce39956 100644
--- a/python/freetoken/layers/quantization/moe/nvfp4.py
+++ b/python/freetoken/layers/quantization/moe/nvfp4.py
@@ -78,6 +78,56 @@ def apply(self, layer, x, topk_weights, topk_ids, view: ExpertView, *, is_prefil
         return fused_experts_decode_nvfp4_marlin(x, *banks, topk_weights, topk_ids, layer.activation, layer.apply_router_weight_on_input, alpha, limit)
 
 
+class W4a4Nvfp4MoEKernel(TritonNvfp4MoEKernel):
+    """The Triton kernels' banks and decode, with the prefill GEMMs on NVFP4 activations (W4A4) where the GPU and the batch favour it.
+
+    Needs a checkpoint whose experts declare static activation scales; the GPU-resident alphas carry them per expert (gate_up_alpha: the larger of gate's and up's input_scale).
+    """
+
+    name = "w4a4"
+
+    def unusable_reason(self, cfg: MoEConfig) -> str | None:
+        reason = super().unusable_reason(cfg)
+        if reason:
+            return reason
+        if cfg.scheme is None or not cfg.scheme.has("input_scale"):
+            return "the experts declare no static activation scale (W4A16 checkpoint)"
+        if cfg.apply_router_weight_on_input:
+            return "the routing weight scales the input, which is quantized once per token"
+        return None
+
+    def worth_it(self, cfg: MoEConfig) -> bool:
+        from freetoken.kernel import nvfp4_moe
+
+        return not backend.is_rocm() and torch.cuda.is_available() and nvfp4_moe.w4a4_available(torch.cuda.current_device())
+
+    def layout(self, cfg: MoEConfig) -> dict[str, BankSpec]:
+        return {
+            **super().layout(cfg),
+            "gate_up_alpha": BankSpec((), torch.float32, resident=True),
+            "down_alpha": BankSpec((), torch.float32, resident=True),
+        }
+
+    def pack(self, pieces, cfg: MoEConfig, out):
+        super().pack(pieces, cfg, out)
+        if not all(f"{p}_input" in pieces for p in ("gate", "up", "down")):
+            raise ValueError("w4a4 nvfp4 MoE kernel: the checkpoint has no expert input_scale tensors")
+        gate, up, down = (pieces[f"{p}_input"].reshape(len(pieces[f"{p}_input"]), -1)[:, 0].float() for p in ("gate", "up", "down"))
+        return {"gate_up_alpha": torch.maximum(gate, up), "down_alpha": down}
+
+    def apply(self, layer, x, topk_weights, topk_ids, view: ExpertView, *, is_prefill: bool):
+        if not is_prefill:
+            return super().apply(layer, x, topk_weights, topk_ids, view, is_prefill=False)
+        from freetoken.moe.fused_nvfp4 import fused_experts_nvfp4_w4a4
+
+        t = view.tensors
+        assert view.alphas is not None
+        return fused_experts_nvfp4_w4a4(
+            x, t["gate_up"], t["gate_up_scale"], t["gate_up_global"], t["down"], t["down_scale"], t["down_global"],
+            view.alphas[0], view.alphas[1], topk_weights, topk_ids, view.n, layer.activation, float(layer.alpha), limit_or_inf(layer),
+        )
+
+
 # ---------------------------------------------------------------------------
 # Marlin: vLLM's fused Marlin MoE (W4A16, dequant in kernel) over Marlin-tiled banks.
 # The Marlin kernels come from vLLM (not sgl-kernel) deliberately: vLLM ships the
@@ -570,4 +620,4 @@ def apply(self, layer, x, topk_weights, topk_ids, view: ExpertView, *, is_prefil
 
 @register_method(QuantKind.NVFP4, LayerKind.MOE)
 class Nvfp4MoEMethod(MoEMethod):
-    candidates = (TritonNvfp4MoEKernel, MarlinNvfp4MoEKernel, B12xNvfp4MoEKernel)
+    candidates = (W4a4Nvfp4MoEKernel, TritonNvfp4MoEKernel, MarlinNvfp4MoEKernel, B12xNvfp4MoEKernel)
diff --git a/python/freetoken/models/nvfp4_banks.py b/python/freetoken/models/nvfp4_banks.py
index 04da05091..e30429c06 100644
--- a/python/freetoken/models/nvfp4_banks.py
+++ b/python/freetoken/models/nvfp4_banks.py
@@ -39,6 +39,28 @@ def _ingest_global(spec: "Nvfp4ExpertSourceSpec", tensor: torch.Tensor) -> torch
     return tensor.to(torch.float16)
 
 
+def _activation_scale():
+    """The checkpoint tensor of the NVFP4 static activation scale under the installed dialect, or None."""
+    from freetoken.layers.quantization import QuantKind, get_quant_config
+
+    try:
+        quant = get_quant_config()
+    except RuntimeError:
+        return None
+    if quant is None or QuantKind.NVFP4 not in quant.STORAGE:
+        return None
+    return quant.stored_tensors(QuantKind.NVFP4).get("input_scale")
+
+
+def _ingest(spec: "Nvfp4ExpertSourceSpec", act, role: str, tensor: torch.Tensor) -> torch.Tensor:
+    if role.endswith("_global"):
+        return _ingest_global(spec, tensor)
+    if role.endswith("_input"):
+        tensor = tensor.float()
+        return 1.0 / tensor if act.reciprocal else tensor
+    return tensor
+
+
 def _num_moe_layers(config) -> int:
     value = getattr(config, "num_moe_layers", None)
     if value is not None:
@@ -77,7 +99,8 @@ def iter_nvfp4_expert_pieces(
 ):
     """One piece per routed expert: ``gate`` / ``up`` / ``down`` codes plus their ``_scale``
     (fp8 block scales) and ``_global`` (the per-tensor scale, reciprocal for quant-side dialects,
-    fp16) companions, straight from the safetensors shards.
+    fp16) companions, straight from the safetensors shards. When every expert also stores the
+    dialect's static activation scale, the pieces carry it as fp32 ``_input`` (dequant side).
 
     Serial reads walk the shards in order; ``parallel`` uses the chunked O_DIRECT reader. Either
     way tensors of one expert may span shards, so they are grouped by (layer, expert) as they land.
@@ -92,6 +115,7 @@ def iter_nvfp4_expert_pieces(
     weight_map = safetensors_weight_map(folder)
 
     keep = None if layers is None else set(layers)
+    act = _activation_scale()
     wanted: dict[str, tuple[int, int, str]] = {}
     for name in weight_map:
         match = spec.key_pattern.match(name)
@@ -107,9 +131,20 @@ def iter_nvfp4_expert_pieces(
         if kind not in ("weight", "weight_scale", "weight_scale_2"):
             raise ValueError(f"{spec.desc}: unknown NVFP4 expert tensor kind {kind!r}")
         wanted[name] = (bank_layer, int(match.group("expert")), spec.proj_to_role[proj] + _kind_suffix(kind))
-    expected = (_num_moe_layers(config) if keep is None else len(keep)) * config.num_experts * 9
-    if len(wanted) != expected:
-        raise ValueError(f"{spec.desc}: found {len(wanted)} expert tensors, expected {expected}")
+        if kind == "weight" and act is not None:
+            act_name = name[: match.start("kind")] + act.name
+            if act_name in weight_map:
+                wanted[act_name] = (bank_layer, int(match.group("expert")), spec.proj_to_role[proj] + "_input")
+    experts = (_num_moe_layers(config) if keep is None else len(keep)) * config.num_experts
+    acts = [name for name, (_, _, role) in wanted.items() if role.endswith("_input")]
+    # a partial set serves no W4A4 kernel; the pieces then carry none
+    if len(acts) != experts * 3:
+        for name in acts:
+            del wanted[name]
+        acts = []
+    per_expert = 12 if acts else 9
+    if len(wanted) != experts * per_expert:
+        raise ValueError(f"{spec.desc}: found {len(wanted)} expert tensors, expected {experts * per_expert}")
 
     def _serial():
         by_shard: dict[str, list[str]] = collections.defaultdict(list)
@@ -121,21 +156,16 @@ def _serial():
             drop(path)
             with safetensors.safe_open(path, framework="pt", device="cpu") as f:
                 for name in by_shard[shard]:
-                    tensor = f.get_tensor(name)
-                    if wanted[name][2].endswith("_global"):
-                        tensor = _ingest_global(spec, tensor)
-                    yield name, tensor
+                    yield name, _ingest(spec, act, wanted[name][2], f.get_tensor(name))
             drop(path)
 
     def _parallel():
         from freetoken.models.weight import iter_expert_tensors_parallel
 
         for name, tensor in iter_expert_tensors_parallel(folder, lambda n: n in wanted, workers=workers, chunk=chunk):
-            if wanted[name][2].endswith("_global"):
-                tensor = _ingest_global(spec, tensor)
-            yield name, tensor
+            yield name, _ingest(spec, act, wanted[name][2], tensor)
 
-    return per_expert_pieces(_parallel() if parallel else _serial(), wanted.get, tensors_per_expert=9)
+    return per_expert_pieces(_parallel() if parallel else _serial(), wanted.get, tensors_per_expert=per_expert)
 
 
 __all__ = ["Nvfp4ExpertSourceSpec", "iter_nvfp4_expert_pieces"]
diff --git a/python/freetoken/moe/expert_pieces.py b/python/freetoken/moe/expert_pieces.py
index 917d84bbc..74c4a41b0 100644
--- a/python/freetoken/moe/expert_pieces.py
+++ b/python/freetoken/moe/expert_pieces.py
@@ -151,7 +151,7 @@ def per_expert_pieces(
     """Group per-expert tensors into one piece per expert, in whatever order they arrive.
 
     ``locate(name)`` -> ``(bank_layer, expert, role)`` or None to skip; a role tensor is one
-    expert's ``[...]`` row and becomes the piece's ``[1, ...]`` entry. Globals may be scalars.
+    expert's ``[...]`` row and becomes the piece's ``[1, ...]`` entry. Globals and activation scales may be scalars.
     """
     pending: dict[tuple[int, int], dict[str, torch.Tensor]] = {}
     for name, tensor in tensors:
@@ -160,7 +160,7 @@ def per_expert_pieces(
             continue
         bank_layer, expert, role = where
         piece = pending.setdefault((bank_layer, expert), {})
-        piece[role] = tensor.reshape(1, -1) if role.endswith("_global") else tensor.unsqueeze(0)
+        piece[role] = tensor.reshape(1, -1) if role.endswith(("_global", "_input")) else tensor.unsqueeze(0)
         if len(piece) == tensors_per_expert:
             del pending[(bank_layer, expert)]
             yield bank_layer, expert, expert + 1, piece
diff --git a/python/freetoken/moe/fused_nvfp4.py b/python/freetoken/moe/fused_nvfp4.py
index 197ef5472..6bedd57b1 100644
--- a/python/freetoken/moe/fused_nvfp4.py
+++ b/python/freetoken/moe/fused_nvfp4.py
@@ -10,6 +10,7 @@
 import torch
 
 from freetoken.kernel import moe_sum_reduce_triton, nvfp4_moe
+from freetoken.kernel.triton.nvfp4_act_quant import quantize_nvfp4
 from freetoken.layers import gated_act_and_mul
 from freetoken.moe.fused import moe_align_block_size
 
@@ -156,8 +157,60 @@ def fused_experts_nvfp4(
     return out
 
 
+def fused_experts_nvfp4_w4a4(
+    hidden_states: torch.Tensor,
+    gate_up_packed: torch.Tensor,
+    gate_up_scale: torch.Tensor,
+    gate_up_global: torch.Tensor,
+    down_packed: torch.Tensor,
+    down_scale: torch.Tensor,
+    down_global: torch.Tensor,
+    gate_up_input_scale: torch.Tensor,
+    down_input_scale: torch.Tensor,
+    topk_weights: torch.Tensor,
+    topk_ids: torch.Tensor,
+    num_experts: int,
+    activation: str = "silu",
+    act_alpha: float = 1.702,
+    act_limit: float = 7.0,
+) -> torch.Tensor:
+    """Prefill NVFP4 MoE with NVFP4 activations (W4A4) where the GPU and the batch favour it, else :func:`fused_experts_nvfp4`.
+
+    ``*_input_scale`` are the per-expert static activation scales of the ``num_experts`` banks; each GEMM input is quantized under their max, as vLLM does.
+    """
+    M, H = hidden_states.shape
+    top_k = topk_ids.shape[1]
+    inter = gate_up_packed.shape[1] // 2
+    cfg = nvfp4_moe.w4a4_prefill_config(M, top_k, num_experts, H, inter, activation, hidden_states.dtype, hidden_states.device.index)
+    if cfg is None:
+        return fused_experts_nvfp4(
+            hidden_states, gate_up_packed, gate_up_scale, gate_up_global, down_packed, down_scale, down_global,
+            topk_weights, topk_ids, num_experts, activation, False, act_alpha, act_limit,
+        )
+    dev = hidden_states.device
+    g = torch.stack((gate_up_input_scale[:num_experts], down_input_scale[:num_experts])).float().amax(dim=1)
+    sorted_ids, expert_ids, ntpp = moe_align_block_size(topk_ids, cfg["BLOCK_SIZE_M"], num_experts)
+    tw = topk_weights.reshape(-1).contiguous()
+    num_valid = topk_ids.numel()
+
+    xq, xs = quantize_nvfp4(hidden_states, g[0:1])
+    q = torch.empty((num_valid, inter // 2), device=dev, dtype=torch.uint8)
+    qs = torch.empty((num_valid, inter // 16), device=dev, dtype=torch.float8_e4m3fn)
+    nvfp4_moe.w4a4_gate_up_gemm(
+        xq, xs, g[0:1], gate_up_packed, gate_up_scale, gate_up_global, q, qs, g[1:2], sorted_ids, expert_ids, ntpp,
+        num_valid, top_k, activation, act_alpha, act_limit, cfg["gate_up"],
+    )
+    out = torch.zeros_like(hidden_states)
+    nvfp4_moe.w4a4_down_gemm(
+        q, qs, g[1:2], down_packed, down_scale, down_global, out, tw, sorted_ids, expert_ids, ntpp, num_valid, top_k,
+        cfg["down"],
+    )
+    return out
+
+
 __all__ = [
     "fused_experts_decode_nvfp4_marlin",
     "fused_experts_decode_nvfp4_serial",
     "fused_experts_nvfp4",
+    "fused_experts_nvfp4_w4a4",
 ]
diff --git a/python/freetoken/moe/legacy_format.py b/python/freetoken/moe/legacy_format.py
index 973a3dba5..8f347e5c3 100644
--- a/python/freetoken/moe/legacy_format.py
+++ b/python/freetoken/moe/legacy_format.py
@@ -30,12 +30,16 @@ def legacy_bank_names(quant_format: str) -> dict[str, str]:
     (QuantKind.NONE, "fused"): "bf16",
     (QuantKind.FP8_BLOCK, "triton"): "fp8_block",
     (QuantKind.NVFP4, "triton"): "nvfp4",
+    # same banks as triton, plus the resident input scales; a "nvfp4" file reads back as triton
+    (QuantKind.NVFP4, "w4a4"): "nvfp4",
     (QuantKind.NVFP4, "marlin"): "nvfp4_marlin",
     (QuantKind.NVFP4, "b12x"): "nvfp4_b12x",
     (QuantKind.MXFP4, "triton_gptoss"): "mxfp4_triton",
     (QuantKind.MXFP4, "triton"): "ds_fp4",
 }
-_KIND_KERNEL = {fmt: kk for kk, fmt in LEGACY_FORMAT.items()}
+_KIND_KERNEL: dict[str, tuple[QuantKind, str]] = {}
+for _kk, _fmt in LEGACY_FORMAT.items():
+    _KIND_KERNEL.setdefault(_fmt, _kk)
 
 
 def legacy_format_for(kind: QuantKind, kernel: str) -> str:
diff --git a/tests/moe/test_nvfp4_backends.py b/tests/moe/test_nvfp4_backends.py
index 29c042ba4..fd98c2ee9 100644
--- a/tests/moe/test_nvfp4_backends.py
+++ b/tests/moe/test_nvfp4_backends.py
@@ -10,13 +10,15 @@
     baseline-kernel equality guard. This is the path used on sm_120 + CUDA 12.x.
   - Marlin (sm_80..sm_99, e.g. H100): prefill + decode + overlap.
   - b12x (sm_120 + CUDA>=13): pack + fused decode, skipped where the kernel cannot run.
-  - Gluon (sm_120): the prefill GEMMs it takes over from the Triton kernel, against that kernel.
+  - Gluon (sm_120): the prefill GEMMs it takes over from the Triton kernel, against that kernel, and the W4A4 prefill
+    against a torch emulation of its NVFP4 activation rounding.
 Banks are built the way the loader builds them: the bound kernel object packs native pieces.
 """
 
 from __future__ import annotations
 
 import importlib.util
+import math
 
 import pytest
 import torch
@@ -413,6 +415,100 @@ def counted(*args):
     torch.testing.assert_close(out, ref, rtol=1e-2, atol=1e-2 * float(ref.abs().max()))
 
 
+def _round_e2m1(v: torch.Tensor) -> torch.Tensor:
+    """Round to the nearest e2m1 value, ties to the even code, saturating at 6 (cvt.rn.satfinite)."""
+    grid = _E2M1[:8].to(v.device)
+    a = v.abs().clamp(max=6.0)
+    hi = torch.searchsorted(grid, a.contiguous()).clamp(max=7)
+    lo = (hi - 1).clamp(min=0)
+    d_lo, d_hi = a - grid[lo], grid[hi] - a
+    code = torch.where(d_lo < d_hi, lo, torch.where(d_hi < d_lo, hi, torch.where(lo % 2 == 0, lo, hi)))
+    return torch.copysign(grid[code], v)
+
+
+def _nvfp4_qdq(x: torch.Tensor, g: float) -> torch.Tensor:
+    """x [R, K] through NVFP4 under the static global g, back to fp32."""
+    xb = x.float().view(x.shape[0], -1, 16)
+    # fp32 scalar math like the kernels: a one-ulp difference in 6 * g can flip an e4m3 block scale
+    g = torch.tensor(g, dtype=torch.float32, device=x.device)
+    sf = ((xb.abs().amax(-1) / (6.0 * g)).to(torch.float8_e4m3fn).float() * g).unsqueeze(-1)
+    inv = torch.where(sf > 0, 1.0 / sf, torch.zeros_like(sf))
+    return (_round_e2m1(xb * inv) * sf).view_as(x)
+
+
+def _gated_ref(h: torch.Tensor, activation: str, inter: int, alpha: float, limit: float) -> torch.Tensor:
+    gate, up = h[:, :inter], h[:, inter:]
+    if activation == "gelu_tanh":
+        return torch.nn.functional.gelu(gate, approximate="tanh") * up
+    if activation in ("swiglu_clamp", "swigluoai"):
+        gate = gate.clamp(max=limit)
+        return gate * torch.sigmoid(alpha * gate) * (up.clamp(-limit, limit) + (activation == "swigluoai"))
+    return torch.nn.functional.silu(gate) * up
+
+
+@cuda
+@pytest.mark.skipif(not torch.cuda.is_available() or torch.cuda.get_device_capability() != (12, 0), reason="needs sm_120")
+@pytest.mark.parametrize("hidden_size,inter,activation,tokens", [
+    (256, 128, "silu", 256), (320, 192, "gelu_tanh", 1024), (512, 1024, "swiglu_clamp", 256), (256, 128, "swigluoai", 256),
+])
+def test_w4a4_prefill_matches_nvfp4_activation_reference(hidden_size, inter, activation, tokens):
+    """The w4a4 kernel from checkpoint pieces to the prefill output, against a torch emulation of the same NVFP4 activation rounding.
+
+    Covers 256- / 128-row gated tiles, K % 256 != 0 in both GEMMs, the persistent down GEMM (I < 1024) and the other one, BLOCK_SIZE_M 64 and 128.
+    """
+    from types import SimpleNamespace
+
+    from freetoken.kernel import nvfp4_moe
+    from freetoken.layers.quantization.moe.base import ExpertView, MoEConfig
+    from freetoken.layers.quantization.moe.nvfp4 import W4a4Nvfp4MoEKernel
+
+    device = torch.device("cuda")
+    if not nvfp4_moe.w4a4_available(device.index or 0):
+        pytest.skip("the W4A4 Gluon kernels need Triton's Gluon on sm_120")
+    alpha, limit = (1.5, 2.0) if activation in ("swiglu_clamp", "swigluoai") else (1.0, None)
+    g = torch.Generator(device=device).manual_seed(0)
+
+    def bank(n, k):
+        packed = torch.randint(0, 256, (E, n, k // 2), dtype=torch.uint8, device=device, generator=g)
+        scale = (torch.rand(E, n, k // 16, device=device, generator=g) + 0.25).to(torch.float8_e4m3fn)
+        return packed, scale, (torch.rand(E, n, device=device, generator=g) * 0.02 + 0.01).to(torch.float16)
+
+    gu, gus, gug = bank(2 * inter, hidden_size)
+    dn, dns, dng = bank(hidden_size, inter)
+    hidden = torch.randn(tokens, hidden_size, dtype=torch.bfloat16, device=device, generator=g)
+    topk_ids = torch.randint(0, E, (tokens, TOPK), dtype=torch.int32, device=device, generator=g)
+    topk_weights = torch.rand(tokens, TOPK, device=device, generator=g)
+
+    # static scales near a calibration's amax / (6 * 448), off by an irrational factor: with amax / 2688 itself many values sit exactly on an
+    # e2m1 midpoint, where the kernel's approximate fp32 reciprocal and torch's exact one round apart; per expert, below the max for some experts
+    g1 = float(hidden.float().abs().max()) / 2688 * math.pi / 3
+    xq = _nvfp4_qdq(hidden, g1)
+    w_gu = torch.stack([_dequant_ref(gu[e], gus[e], gug[e]) for e in range(E)])
+    w_dn = torch.stack([_dequant_ref(dn[e], dns[e], dng[e]) for e in range(E)])
+    routes = topk_ids.reshape(-1).long()
+    h = torch.bmm(w_gu[routes], xq.repeat_interleave(TOPK, 0).unsqueeze(-1)).squeeze(-1).bfloat16().float()
+    act = _gated_ref(h, activation, inter, alpha, limit or 0.0).bfloat16().float()
+    g2 = float(act.abs().max()) / 2688 * math.pi / 3
+    ref = torch.bmm(w_dn[routes], _nvfp4_qdq(act, g2).unsqueeze(-1)).squeeze(-1)
+    ref = (ref * topk_weights.reshape(-1, 1)).view(tokens, TOPK, hidden_size).sum(1)
+
+    def input_scale(top):
+        return torch.full((E, 1), top, device=device).index_fill_(0, torch.tensor([1, 2], device=device), top / 2)
+
+    pieces = {"gate_up": gu, "gate_up_scale": gus, "gate_up_global": gug, "down": dn, "down_scale": dns, "down_global": dng,
+              "gate_input": input_scale(g1), "up_input": input_scale(g1).roll(1, 0), "down_input": input_scale(g2)}
+    cfg = MoEConfig(num_experts=E, hidden=hidden_size, intermediate=inter, top_k=TOPK, activation=activation, alpha=alpha, limit=limit)
+    kernel = W4a4Nvfp4MoEKernel()
+    layout = kernel.layout(cfg)
+    banks = {role: torch.empty((E, *spec.shape), dtype=spec.dtype, device=device) for role, spec in layout.items() if not spec.resident}
+    alphas = kernel.pack(pieces, cfg, banks)
+    view = ExpertView(banks, n=E, alphas=(alphas["gate_up_alpha"], alphas["down_alpha"]))
+    layer = SimpleNamespace(activation=activation, alpha=alpha, limit=limit, apply_router_weight_on_input=False)
+    out = kernel.apply(layer, hidden, topk_weights, topk_ids, view, is_prefill=True)
+    assert torch.equal(alphas["gate_up_alpha"], torch.maximum(pieces["gate_input"], pieces["up_input"]).squeeze(1))
+    torch.testing.assert_close(out.float(), ref, rtol=1e-2, atol=1e-2 * float(ref.abs().max()))
+
+
 def _triton_cache(device, *, cache_size=S, prefill_overlap=False):
     """Native 6-bank NVFP4 cache (no repack), consumed directly by the Triton kernels.
     The banks are not transformed, so ``sources`` doubles as the dequant reference."""