diff --git a/docs/kernel-contract.md b/docs/kernel-contract.md
index 0d8241d5c2..2dd2a3a5e4 100644
--- a/docs/kernel-contract.md
+++ b/docs/kernel-contract.md
@@ -171,7 +171,8 @@ RADIO's fused bias, activation and residual (`tf_vit_bias_act`) replaces `tf_vit
 CUDA also checks GLM-5.3-Flash's KV page kernels on a page table out of order: the latent write and the pooled keys
 (`tf_glm_lat_write`, `tf_glm_pool_keys`), the index scores and the selected tokens' rows (`tf_glm_scores`,
 `tf_glm_dsa_tokens`), and the MLA chunks through their merge (`tf_glm_chunks_16` and `tf_glm_chunks_32`, dense and
-sparse).
+sparse). Its prompt indexer (`tf_glm_index_scores`, `tf_glm_index_pools`) is bit-equal to `tf_glm_scores` and
+`tf_glm_select_rows`; the pools' stated rule ranks any NaN above +inf, as the GPU's `s + 0` makes every NaN one positive NaN.
 
 GLM's Metal window projections register narrow and wide rows at 2, 3, 4, 5, 6 and 8 bits, plus both gathered paths.
 Independent FP64 linear cases accompany byte replacements against frozen pre-window-optimization kernels.
@@ -191,8 +192,8 @@ and its split-K form, prompt16's FP8 tile, qmmf's FP4 32-row cluster tile, five
 and FP32-output paths, FP8G's 32-row cluster tile with FP32 output, the common EXL3 prompt products, Super's
 decode-sized latent top-k and NVFP4 up experts, the Mamba scan and Super's NVFP4 prompt down kernels, three of
 GLM-5.3-Flash's packed matmuls and three of its DFlash2 launches, and the training launches it compiles differently,
-and nvcc 13.3 keeps 13.2's and changes `tf_nemo_topk` on sm_120 and sm_121 and four of GLM-5.3-Flash's KV page
-launches. `kernel-contract` and
+and nvcc 13.3 keeps 13.2's and changes `tf_nemo_topk` on sm_120 and sm_121, four of GLM-5.3-Flash's KV page
+launches and its two prompt-indexer launches. `kernel-contract` and
 `kernel-contract-limits` read the release from each image's `.note.nv.tkinfo`, check the image against the newest set
 at or below that release (13.0's when an image names none), and print which set they used. A release with no set of
 its own is checked against the one below it, so its differences show as stale records: re-derive them with
diff --git a/tools/glm_index_emulate/cuda_bf16.h b/tools/glm_index_emulate/cuda_bf16.h
new file mode 100644
index 0000000000..61672c3aa6
--- /dev/null
+++ b/tools/glm_index_emulate/cuda_bf16.h
@@ -0,0 +1,16 @@
+// bfloat16 for the host emulation: the type and its two conversions.
+#pragma once
+#include <cstdint>
+#include <cstring>
+struct __nv_bfloat16 { uint16_t x; };
+inline float __bfloat162float(__nv_bfloat16 v) {
+  const uint32_t b = (uint32_t)v.x << 16;
+  float f;
+  std::memcpy(&f, &b, 4);
+  return f;
+}
+inline __nv_bfloat16 __float2bfloat16_rn(float f) {
+  uint32_t b;
+  std::memcpy(&b, &f, 4);
+  return {(uint16_t)((b + 0x7FFF + ((b >> 16) & 1)) >> 16)};
+}
diff --git a/tools/glm_index_emulate/cuda_runtime.h b/tools/glm_index_emulate/cuda_runtime.h
new file mode 100644
index 0000000000..918ab0b427
--- /dev/null
+++ b/tools/glm_index_emulate/cuda_runtime.h
@@ -0,0 +1,62 @@
+// Host stand-ins for the CUDA names glm_index_*.cuh use: a block's threads are OS threads, __syncthreads a barrier.
+#pragma once
+#include <algorithm>
+#include <atomic>
+#include <barrier>
+#include <cmath>
+#include <cstdint>
+#include <cstring>
+#define __global__
+#define __device__
+#define __forceinline__ inline
+#define __restrict__
+#define __launch_bounds__(x)
+#define __shared__ static
+#define __align__(x) __attribute__((aligned(x)))
+struct dim3 { unsigned x, y, z; };
+struct float4 { float x, y, z, w; };
+inline thread_local dim3 threadIdx, blockIdx;
+namespace emu {
+inline std::barrier<>* block_bar;
+inline std::barrier<>* warp_bar[32];
+inline uint32_t xbuf[1024];
+inline const void* pbuf[1024];
+}  // namespace emu
+inline void __syncthreads() { emu::block_bar->arrive_and_wait(); }
+inline void warp_sync() { emu::warp_bar[threadIdx.x >> 5]->arrive_and_wait(); }
+template <class T>
+inline T shuffle(T v, int src_of(int lane)) {
+  const int t = (int)threadIdx.x, l = t & 31, base = t - l;
+  uint32_t w;
+  std::memcpy(&w, &v, 4);
+  emu::xbuf[t] = w;
+  warp_sync();
+  const int s = src_of(l);
+  std::memcpy(&v, &emu::xbuf[base + (s >= 0 ? s : l)], 4);
+  warp_sync();
+  return v;
+}
+inline int __shfl_up_sync(unsigned, int v, int d) {
+  static thread_local int delta;
+  delta = d;
+  return shuffle(v, [](int l) { return l >= delta ? l - delta : -1; });
+}
+inline float __shfl_xor_sync(unsigned, float v, int o) {
+  static thread_local int mask;
+  mask = o;
+  return shuffle(v, [](int l) { return l ^ mask; });
+}
+inline int atomicAdd(int* p, int v) { return std::atomic_ref<int>(*p).fetch_add(v); }
+inline int atomicMax(int* p, int v) {
+  std::atomic_ref<int> a(*p);
+  int old = a.load();
+  while (old < v && !a.compare_exchange_weak(old, v)) {}
+  return old;
+}
+inline int __float_as_int(float f) { int i; std::memcpy(&i, &f, 4); return i; }
+inline float __fadd_rn(float a, float b) { return a + b; }
+inline float __fmul_rn(float a, float b) { return a * b; }
+inline float __fmaf_rn(float a, float b, float c) { return std::fma(a, b, c); }
+using std::fmaxf;
+using std::max;
+using std::min;
diff --git a/tools/glm_index_emulate/emulate.cpp b/tools/glm_index_emulate/emulate.cpp
new file mode 100644
index 0000000000..7a414b9c91
--- /dev/null
+++ b/tools/glm_index_emulate/emulate.cpp
@@ -0,0 +1,225 @@
+// glm_index_*.cuh on the host, a block's threads at once, against plain models (zig build test-glm-index-emulate)
+#define GLM_INDEX_HOST
+#include <cstdio>
+#include <functional>
+#include <random>
+#include <thread>
+#include <vector>
+#include "cuda_bf16.h"
+#include "cuda_runtime.h"
+
+namespace glm_index {
+inline void cp16(void* d, const void* s, bool ok) { ok ? (void)std::memcpy(d, s, 16) : (void)std::memset(d, 0, 16); }
+inline void cp_commit() {}
+inline void cp_wait() {}
+inline void ldm4(uint32_t* r, const void* p) {  // matrix i, row j of the 8 x 8 b16 tiles at the address of lane 8 i + j
+  const int t = (int)threadIdx.x, l = t & 31, base = t - l;
+  emu::pbuf[t] = p;
+  warp_sync();
+  for (int i = 0; i < 4; ++i) std::memcpy(&r[i], (const char*)emu::pbuf[base + 8 * i + l / 4] + 4 * (l % 4), 4);
+  warp_sync();
+}
+inline float bf_at(uint32_t w, int half) { return __bfloat162float({(uint16_t)(half ? w >> 16 : w & 0xFFFF)}); }
+inline void mma16816(float* c, const uint32_t* a, const uint32_t* b) {  // m16n8k16, c += a b as one double sum
+  static uint32_t xa[1024][6];
+  const int t = (int)threadIdx.x, l = t & 31, base = t - l, g = l / 4, q = l % 4;
+  for (int i = 0; i < 4; ++i) xa[t][i] = a[i];
+  xa[t][4] = b[0];
+  xa[t][5] = b[1];
+  warp_sync();
+  float out[4];
+  for (int j = 0; j < 4; ++j) {
+    const int row = g + 8 * (j / 2), col = 2 * q + j % 2;
+    double sum = 0;
+    for (int k = 0; k < 16; ++k) {
+      const uint32_t aw = xa[base + (row % 8) * 4 + (k % 8) / 2][(row < 8 ? 0 : 1) + (k < 8 ? 0 : 2)];
+      const uint32_t bw = xa[base + col * 4 + (k % 8) / 2][4 + (k < 8 ? 0 : 1)];
+      sum += (double)bf_at(aw, k % 2) * (double)bf_at(bw, k % 2);
+    }
+    out[j] = (float)((double)c[j] + sum);
+  }
+  warp_sync();
+  for (int j = 0; j < 4; ++j) c[j] = out[j];
+}
+}  // namespace glm_index
+
+#include "../../zig/kernels/cuda/glm_index_pools.cuh"
+#include "../../zig/kernels/cuda/glm_index_scores.cuh"
+
+// every block of the grid in turn, its threads at once
+static void run_grid(unsigned gx, unsigned gy, int threads, const std::function<void()>& kernel) {
+  for (unsigned by = 0; by < gy; by++)
+    for (unsigned bx = 0; bx < gx; bx++) {
+      std::barrier<> block(threads);
+      std::vector<std::barrier<>*> warps;
+      for (int w = 0; w < threads / 32; w++) warps.push_back(new std::barrier<>(32));
+      emu::block_bar = &block;
+      for (int w = 0; w < threads / 32; w++) emu::warp_bar[w] = warps[w];
+      std::vector<std::thread> ts;
+      for (int t = 0; t < threads; t++)
+        ts.emplace_back([=, &kernel] {
+          blockIdx = {bx, by, 0};
+          threadIdx = {(unsigned)t, 0, 0};
+          kernel();
+        });
+      for (auto& t : ts) t.join();
+      for (auto* w : warps) delete w;
+    }
+}
+
+static int fails = 0;
+static uint32_t bits(float f) { uint32_t u; std::memcpy(&u, &f, 4); return u; }
+static float bf(__nv_bfloat16 v) { return __bfloat162float(v); }
+
+struct Case { int64_t pos, rows, np; int loop, heads; bool paged; };
+
+// tf_glm_scores' result for one (row, pool), the model of the mma: a chain of double sums per k16 step, then its rules
+static float model_score(const std::vector<__nv_bfloat16>& q, const std::vector<__nv_bfloat16>& w, int64_t w_stride,
+                         const __nv_bfloat16* key, int64_t r, int heads, float scale, float wscale) {
+  float x[32] = {}, wv[32] = {};
+  for (int h = 0; h < heads; h++) {
+    float acc = 0;
+    for (int kq = 0; kq < 8; kq++) {
+      double sum = 0;
+      for (int k = 0; k < 16; k++) sum += (double)bf(q[(r * heads + h) * 128 + kq * 16 + k]) * (double)bf(key[kq * 16 + k]);
+      acc = (float)((double)acc + sum);
+    }
+    x[h] = std::fmax(acc * scale, 0.0f);
+    wv[h] = bf(w[r * w_stride + h]) * wscale;
+  }
+  float s[8];
+  for (int g = 0; g < 8; g++) {
+    float v = wv[g + 8] * x[g + 8];
+    v = std::fma(wv[g], x[g], v);
+    v = std::fma(wv[g + 16], x[g + 16], v);
+    s[g] = std::fma(wv[g + 24], x[g + 24], v);
+  }
+  float t[8], u[8];
+  for (int g = 0; g < 8; g++) t[g] = s[g] + s[g ^ 4];
+  for (int g = 0; g < 8; g++) u[g] = t[g] + t[g ^ 2];
+  return u[0] + u[1];
+}
+
+static void scores_case(std::mt19937_64& g, const Case& c) {
+  const int D = 128, W = 40;  // head weights in a row of W columns
+  const int64_t visible = (c.pos + c.rows) / 4;
+  const int64_t nb = (std::min(c.np, std::max<int64_t>(512, visible)) + 63) / 64;
+  std::normal_distribution<float> nd(0.0f, 1.0f);
+  std::vector<__nv_bfloat16> q(c.rows * c.heads * D), w(c.rows * W);
+  for (auto& x : q) x = __float2bfloat16_rn(nd(g));
+  for (auto& x : w) x = __float2bfloat16_rn(nd(g));
+  const int64_t pool_rows = (c.pos + c.rows) / 4 + 1;
+  std::vector<__nv_bfloat16> keys(std::max<int64_t>(c.np, pool_rows + 200) * D);
+  for (auto& x : keys) x = __float2bfloat16_rn(nd(g) * 2);
+  // pages of 128 pools, shuffled, when paged: the table maps a logical page to its place
+  const int64_t npages = (keys.size() / D + 127) / 128, page_elems = 128 * D;
+  std::vector<int> pt(npages), order(npages);
+  for (int i = 0; i < npages; i++) order[i] = i;
+  std::shuffle(order.begin(), order.end(), g);
+  std::vector<__nv_bfloat16> store(npages * page_elems);
+  for (int64_t p = 0; p < (int64_t)(keys.size() / D); p++) {
+    pt[p / 128] = order[p / 128];
+    std::memcpy(&store[(int64_t)order[p / 128] * page_elems + (p % 128) * D], &keys[p * D], D * 2);
+  }
+  const __nv_bfloat16* pk = c.paged ? store.data() : keys.data();
+  const float scale = 0.0883883461f, wscale = 0.176776692f;
+  const int pos = (int)c.pos;
+  std::vector<float> out(c.rows * c.np);
+  std::vector<uint32_t> raw(out.size(), 0x7FC0DEAD);
+  std::memcpy(out.data(), raw.data(), out.size() * 4);
+  run_grid((c.rows + 7) / 8, (nb + c.loop - 1) / c.loop, 256, [&] {
+    tf_glm_index_scores(q.data(), w.data(), W, pk, out.data(), &pos, c.rows, c.np, nb, scale, wscale, c.heads, c.loop,
+                        c.paged ? pt.data() : nullptr, c.paged ? page_elems * 2 : 0);
+  });
+  int64_t bad = 0, unread = 0;
+  for (int64_t r = 0; r < c.rows; r++)
+    for (int64_t p = 0; p < c.np; p++) {
+      const uint32_t have = bits(out[r * c.np + p]);
+      if (p >= nb * 64) {
+        unread += have != 0x7FC0DEAD;  // never written
+        continue;
+      }
+      const int64_t npool = (c.pos + r + 1) / 4;
+      const float want = p < npool ? model_score(q, w, W, &keys[p * D], r, c.heads, scale, wscale) : -INFINITY;
+      bad += have != bits(want);
+    }
+  printf("scores pos %lld rows %lld np %lld loop %d heads %d %s: %lld differ, %lld untouched columns written\n",
+         (long long)c.pos, (long long)c.rows, (long long)c.np, c.loop, c.heads, c.paged ? "paged" : "plain",
+         (long long)bad, (long long)unread);
+  fails += bad != 0 || unread != 0;
+}
+
+static int tile_cap(int np) { return np <= 12288 ? 1024 : np <= 32768 ? 2048 : 4096; }
+static int tile_target(int np) { return np <= 12288 ? 724 : np <= 32768 ? 1024 : 2048; }
+static uint32_t key(float s) { return glm_index::order_key(s); }
+
+static std::vector<int64_t> reference(const float* row, int n, int k) {
+  std::vector<int> o(n);
+  for (int i = 0; i < n; i++) o[i] = i;
+  std::stable_sort(o.begin(), o.end(), [&](int a, int b) { return key(row[a]) > key(row[b]); });
+  std::vector<int64_t> best(o.begin(), o.begin() + k);
+  std::sort(best.begin(), best.end());
+  return best;
+}
+
+enum Kind { normal, ties, equal, zeros, skewed, decoy };
+
+static void fill(std::mt19937_64& g, float* row, int np, int limit, Kind kind) {
+  std::normal_distribution<float> nd(0.0f, 3.0f);
+  const int vis = limit;  // the row's visible pools; the rest read as stale 1e9, which the kernel must never touch
+  for (int i = 0; i < np; i++) {
+    float v = nd(g);
+    if (kind == ties) v = (float)(int)(g() % 9) * 0.5f;
+    if (kind == equal) v = 1.25f;
+    if (kind == zeros) v = g() % 3 == 0 ? 0.0f : g() % 2 ? -0.0f : (float)(g() % 1000) * 1e-3f - 0.6f;
+    if (kind == skewed) v = i % 997 == 0 ? 100.0f + (float)(g() % 5) : -INFINITY;
+    if (kind == decoy) v = (float)(g() % 1000) * 1e-3f;
+    row[i] = i < vis ? v : 1e9f;
+  }
+  if (kind == decoy && limit > 4096)  // the floor's sample runs high: fewer than K pools at or above it
+    for (int s = 0; s < glm_index::NS * 8; s++) {
+      const int64_t j = (int64_t)(s / 8) * limit / glm_index::NS / 8 * 8 + s % 8;
+      if (j < limit) row[j] = 10.0f + (float)(g() % 100000) * 1e-4f;
+    }
+}
+
+int main() {
+  std::mt19937_64 g(0x101);
+  for (const Case& c : {Case{0, 1, 1024, 2, 32, false}, Case{1000, 19, 1024, 4, 32, false}, Case{3000, 24, 1024, 16, 32, true},
+                        Case{3900, 17, 1000, 3, 32, false}, Case{5000, 9, 2048, 5, 20, true},
+                        Case{6000, 24, 2048, 8, 32, false}, Case{12000, 8, 4096, 3, 32, true}})
+    scores_case(g, c);
+
+  const int K = 512, rows = 2;
+  int paths[3] = {0, 0, 0};
+  for (int np : {1024, 4097, 12289, 32769, 70001})
+    for (int mode = 0; mode < 3; mode++)
+      for (Kind kind : {normal, ties, equal, zeros, skewed, decoy}) {
+        // the first row's position: few visible pools (limit K), about half, or nearly all
+        const int pos = mode == 0 ? 300 : mode == 1 ? 2 * np : 4 * np - rows - 5;
+        std::vector<float> s((size_t)rows * np + 4);
+        float* base = s.data() + (np % 4 ? 1 : 0);  // an unaligned row now and then: the scalar loads
+        std::vector<int> limit(rows);
+        for (int r = 0; r < rows; r++) {
+          limit[r] = (int)std::min<int64_t>(np, std::max<int64_t>((pos + r + 1) / 4, K));
+          fill(g, base + (size_t)r * np, np, limit[r], kind);
+        }
+        std::vector<int64_t> out((size_t)rows * K, -7);
+        std::vector<int> path(rows, -1);
+        run_grid(rows, 1, glm_index::T, [&] {
+          tf_glm_index_pools(base, out.data(), np, &pos, K, tile_cap(np), tile_target(np), path.data());
+        });
+        for (int r = 0; r < rows; r++) {
+          paths[path[r]]++;
+          const auto want = reference(base + (size_t)r * np, limit[r], K);
+          if (!std::equal(want.begin(), want.end(), out.begin() + (size_t)r * K)) {
+            fails++;
+            printf("pools np %d mode %d kind %d row %d (path %d) differs\n", np, mode, (int)kind, r, path[r]);
+          }
+        }
+      }
+  printf("pools paths: whole row %d, candidates %d, whole row after candidates %d\n", paths[0], paths[1], paths[2]);
+  if (paths[1] == 0 || paths[2] == 0) { fails++; printf("a path was never taken\n"); }
+  printf("%s (%d failures)\n", fails ? "FAIL" : "PASS", fails);
+  return fails != 0;
+}
diff --git a/zig/build/cuda_kernels.zig b/zig/build/cuda_kernels.zig
index a84e9bbb9b..887b2d84ba 100644
--- a/zig/build/cuda_kernels.zig
+++ b/zig/build/cuda_kernels.zig
@@ -101,6 +101,7 @@ pub const kernels = [_]Kernel{
     .{ .name = "glm_kda", .flags = &.{ "-O3", "--fmad=false" } }, // GLM: the Python family's glm5_next/cuda/kda.cu
     .{ .name = "glm_glue", .flags = &.{"-O3"} }, // ours: the Python family's Triton glue by hand, its PTX op for op
     .{ .name = "glm_ops", .flags = &.{"-O3"} }, // ours: GLM's integer work the Python family did in torch (the DSA token lists)
+    .{ .name = "glm_index", .flags = &.{"-O3"} }, // a prompt chunk's DSA scores and top-512 pools
     .{ .name = "glm_quant", .flags = torch_ops }, // ours: GLM's draft-only weight quantizers, torch at load in the Python family
     .{ .name = "exl3_linear", .flags = &.{ "-O3", "--expt-relaxed-constexpr" } }, // exl3/linear.cu's device code, tensorfold_exl3_linear_v5
     .{ .name = "exl3_experts", .flags = &.{ "-O3", "--expt-relaxed-constexpr" } }, // exl3/experts.cu's device code, tensorfold_exl3_experts
diff --git a/zig/build/glm_family_host.zig b/zig/build/glm_family_host.zig
index 105ba82eff..766db3392f 100644
--- a/zig/build/glm_family_host.zig
+++ b/zig/build/glm_family_host.zig
@@ -8,6 +8,16 @@ pub fn steps(b: *std.Build, cuda: *std.Build.Module, glm: *std.Build.Module, mod
     b.step("test-glm-startup-host", "Page refusals propagate through rank liveness before engine construction").dependOn(&startup_run.step);
     step.dependOn(&startup_run.step);
     b.step("test-glm-tp3-host", "Uneven KDA head shares through actual projection and load planning, without devices").dependOn(&b.addRunArtifact(b.addTest(.{ .root_module = glm, .filters = &.{"GLM TP3:"} })).step);
+    // TF_GLM_INDEX_CUDA: the launch rules, and (its own slower step) glm_index_*.cuh as host C++ against plain models
+    const index_run = b.addRunArtifact(b.addTest(.{ .root_module = glm, .filters = &.{"GLM index:"} }));
+    b.step("test-glm-index-host", "GLM's prompt indexer launch rules without devices").dependOn(&index_run.step);
+    step.dependOn(&index_run.step);
+    const index_emulate = b.addExecutable(.{ .name = "glm-index-emulate", .root_module = b.createModule(.{ .target = host, .optimize = .ReleaseSafe, .link_libcpp = true }) });
+    index_emulate.root_module.addCSourceFile(.{ .file = b.path("tools/glm_index_emulate/emulate.cpp"), .flags = &.{ "-std=c++20", "-O1" } });
+    index_emulate.root_module.addIncludePath(b.path("tools/glm_index_emulate"));
+    const index_emulated = b.addRunArtifact(index_emulate);
+    index_emulated.expectExitCode(0);
+    b.step("test-glm-index-emulate", "GLM's prompt indexer kernels as host C++ against a plain top-K and a model of tf_glm_scores").dependOn(&index_emulated.step);
     const kv8_run = b.addRunArtifact(b.addTest(.{ .root_module = glm, .filters = &.{"GLM FP8 KV:"} }));
     b.step("test-glm-fp8-kv-host", "GLM latent cache layout, mode agreement and recorded CUDA launches without devices").dependOn(&kv8_run.step);
     step.dependOn(&kv8_run.step);
diff --git a/zig/kernels/cuda/glm_index.cu b/zig/kernels/cuda/glm_index.cu
new file mode 100644
index 0000000000..128f942821
--- /dev/null
+++ b/zig/kernels/cuda/glm_index.cu
@@ -0,0 +1,3 @@
+// GLM-5.3-Flash's prompt-chunk DSA indexer in CUDA (cuda_index.zig): pool scores, then each row's top-K pools.
+#include "glm_index_scores.cuh"
+#include "glm_index_pools.cuh"
diff --git a/zig/kernels/cuda/glm_index_pools.cuh b/zig/kernels/cuda/glm_index_pools.cuh
new file mode 100644
index 0000000000..4882b760ce
--- /dev/null
+++ b/zig/kernels/cuda/glm_index_pools.cuh
@@ -0,0 +1,218 @@
+// GLM-5.3-Flash DSA indexer, a prompt chunk's top-K pools in one pass: a sampled floor, the survivors, a radix select.
+#pragma once
+
+#include <cuda_runtime.h>
+#include <stdint.h>
+
+namespace glm_index {
+
+constexpr int T = 256, WARPS = T / 32, PER = 8, CHUNK = T * PER;  // threads, warps, scores a thread a step, a step
+constexpr int CAP_MAX = 4096, NS = 512;                           // candidates in shared memory, sample runs of 8
+
+struct Smem {
+  uint32_t key[CAP_MAX];  // candidates' order keys, in pool order
+  int idx[CAP_MAX];       // their pools
+  int hist[256];
+  int warp[2][WARPS];  // block_scan's warp totals, two buffers taken in turn
+  int digit, need;     // a radix step's results
+};
+
+// sparse._order_key: a score's bits as a uint32 whose unsigned order is the scores' (-0 counted as +0).
+__device__ __forceinline__ uint32_t order_key(float s) {
+  const int bits = __float_as_int(__fadd_rn(s, 0.0f));
+  return (uint32_t)(bits ^ ((bits >> 31) | (int)0x80000000u));
+}
+
+// Exclusive sum of v over the block's threads in thread order; `total` the block's sum. Every thread calls it.
+__device__ __forceinline__ int block_scan(Smem& sm, int v, int& phase, int& total) {
+  const int lane = threadIdx.x & 31, w = threadIdx.x >> 5, buf = phase & 1;
+  phase++;
+  int x = v;
+#pragma unroll
+  for (int d = 1; d < 32; d <<= 1) {
+    const int y = __shfl_up_sync(0xffffffffu, x, d);
+    if (lane >= d) x += y;
+  }
+  if (lane == 31) sm.warp[buf][w] = x;
+  __syncthreads();
+  int before = 0;
+  total = 0;
+#pragma unroll
+  for (int i = 0; i < WARPS; ++i) {
+    const int s = sm.warp[buf][i];
+    if (i < w) before += s;
+    total += s;
+  }
+  return before + x - v;
+}
+
+// The row's scores: element i's key and pool, and the keys of elements i0 .. i0 + PER - 1 (0 past n).
+struct Row {
+  const float* s;
+  int n;
+  __device__ __forceinline__ uint32_t key(int i) const { return order_key(s[i]); }
+  __device__ __forceinline__ int pool(int i) const { return i; }
+  __device__ __forceinline__ void keys(int i0, uint32_t (&u)[PER]) const {
+    if (i0 + PER <= n && ((uintptr_t)(s + i0) & 15) == 0) {
+      const float4 a = *reinterpret_cast<const float4*>(s + i0), b = *reinterpret_cast<const float4*>(s + i0 + 4);
+      const float v[PER] = {a.x, a.y, a.z, a.w, b.x, b.y, b.z, b.w};
+#pragma unroll
+      for (int j = 0; j < PER; ++j) u[j] = order_key(v[j]);
+    } else {
+#pragma unroll
+      for (int j = 0; j < PER; ++j) u[j] = i0 + j < n ? order_key(s[i0 + j]) : 0u;
+    }
+  }
+};
+
+// The candidates in shared memory, in the same shape.
+struct Cand {
+  const Smem* sm;
+  int n;
+  __device__ __forceinline__ uint32_t key(int i) const { return sm->key[i]; }
+  __device__ __forceinline__ int pool(int i) const { return sm->idx[i]; }
+  __device__ __forceinline__ void keys(int i0, uint32_t (&u)[PER]) const {
+#pragma unroll
+    for (int j = 0; j < PER; ++j) u[j] = i0 + j < n ? sm->key[i0 + j] : 0u;
+  }
+};
+
+// sparse._radix_step from sm.hist: the pass's digit of the need-th best key under prefix, and need less those above.
+__device__ __forceinline__ void radix_step(Smem& sm, int& need, uint32_t& prefix, uint32_t& fixed, int p, int& phase) {
+  const int b = T - 1 - threadIdx.x;  // thread t holds bin 255 - t
+  const int h = sm.hist[b];
+  int total;
+  const int above = block_scan(sm, h, phase, total);  // pools in bins above b
+  if (above + h >= need) atomicMax(&sm.digit, b);
+  __syncthreads();
+  const int hit = sm.digit;  // -1: no bin reaches need, as tf_glm_select_rows takes digit 0 and need less every pool
+  if (hit >= 0 && b == hit) sm.need = need - above;
+  __syncthreads();
+  const uint32_t digit = hit >= 0 ? (uint32_t)hit : 0u;
+  need = hit >= 0 ? sm.need : need - total;
+  prefix |= digit << (24 - 8 * p);
+  fixed |= 0xFFu << (24 - 8 * p);
+}
+
+// A pass's histogram cleared and the digit reset; the caller's previous reads of them are done first.
+__device__ __forceinline__ void hist_clear(Smem& sm) {
+  __syncthreads();
+  sm.hist[threadIdx.x] = 0;
+  if (threadIdx.x == 0) sm.digit = -1;
+  __syncthreads();
+}
+
+// sparse._select_rows over src's n elements: the K best's pools ascending into out, ties to the lower pool.
+template <class Src>
+__device__ void select(Smem& sm, const Src& src, int K, int64_t* __restrict__ out, int& phase) {
+  uint32_t prefix = 0, fixed = 0;
+  int need = K;
+  for (int p = 0; p < 4; ++p) {
+    hist_clear(sm);
+    const int sh = 24 - 8 * p;
+    for (int i = threadIdx.x; i < src.n; i += T) {
+      const uint32_t u = src.key(i);
+      if ((u & fixed) == prefix) atomicAdd(&sm.hist[(u >> sh) & 0xFFu], 1);
+    }
+    __syncthreads();
+    radix_step(sm, need, prefix, fixed, p, phase);
+  }
+  // in pool order: every key above the K-th best's, and the first `need` equal to it, each at its rank among those
+  int gt_run = 0, eq_run = 0;
+  for (int c0 = 0; c0 < src.n; c0 += CHUNK) {
+    const int i0 = c0 + (int)threadIdx.x * PER;
+    uint32_t u[PER];
+    src.keys(i0, u);
+    int g = 0, e = 0;
+#pragma unroll
+    for (int j = 0; j < PER; ++j) {
+      const bool ok = i0 + j < src.n;
+      g += ok && u[j] > prefix;
+      e += ok && u[j] == prefix;
+    }
+    int total;
+    const int ex = block_scan(sm, (g << 16) | e, phase, total);
+    int gb = gt_run + (ex >> 16), eb = eq_run + (ex & 0xFFFF);
+#pragma unroll
+    for (int j = 0; j < PER; ++j) {
+      if (i0 + j >= src.n) break;
+      if (u[j] > prefix) {
+        out[gb + min(eb, need)] = (int64_t)src.pool(i0 + j);
+        gb++;
+      } else if (u[j] == prefix) {
+        if (eb < need) out[gb + eb] = (int64_t)src.pool(i0 + j);
+        eb++;
+      }
+    }
+    gt_run += total >> 16;
+    eq_run += total & 0xFFFF;
+  }
+}
+
+// The 16-bit key prefix of the rank of NS runs of 8 scores at even spacing that scales to TARGET pools. Speed only.
+__device__ uint32_t sample_floor(Smem& sm, const float* row, int NP, int TARGET, int& phase) {
+  uint32_t prefix = 0, fixed = 0;
+  int need = max((int)((int64_t)TARGET * (NS * 8) / NP), 1);
+  for (int p = 0; p < 2; ++p) {
+    hist_clear(sm);
+    const int sh = 24 - 8 * p;
+    for (int s = threadIdx.x; s < NS * 8; s += T) {
+      const int64_t j = (int64_t)(s / 8) * NP / NS / 8 * 8 + s % 8;
+      const uint32_t u = order_key(j < NP ? row[j] : 0.0f);
+      if ((u & fixed) == prefix) atomicAdd(&sm.hist[(u >> sh) & 0xFFu], 1);
+    }
+    __syncthreads();
+    radix_step(sm, need, prefix, fixed, p, phase);
+  }
+  return prefix;
+}
+
+}  // namespace glm_index
+
+// tf_glm_select_rows' pools for scores [R][np]; CAP (<= 4096), TARGET speed only; PATH null or [R] int32; grid (R).
+extern "C" __global__ void __launch_bounds__(glm_index::T)
+    tf_glm_index_pools(const float* __restrict__ S, int64_t* __restrict__ OUT, int64_t np, const int* __restrict__ pos,
+                       int K, int CAP, int TARGET, int* __restrict__ PATH) {
+  using namespace glm_index;
+  __shared__ Smem sm;
+  CAP = min(CAP, CAP_MAX);
+  const int64_t r = blockIdx.x, vis = ((int64_t)pos[0] + r + 1) / 4;
+  const float* row = S + r * np;
+  int64_t* out = OUT + r * K;
+  const int NP = (int)min(np, max(vis, (int64_t)K));
+  int phase = 0;
+  if (NP <= CAP) {
+    select(sm, Row{row, NP}, K, out, phase);
+    if (PATH != nullptr && threadIdx.x == 0) PATH[r] = 0;
+    return;
+  }
+  const uint32_t lo = sample_floor(sm, row, NP, TARGET, phase);
+  const Row src{row, NP};
+  int n = 0;
+  for (int c0 = 0; c0 < NP; c0 += CHUNK) {
+    const int i0 = c0 + (int)threadIdx.x * PER;
+    uint32_t u[PER];
+    src.keys(i0, u);
+    int m = 0;
+#pragma unroll
+    for (int j = 0; j < PER; ++j) m += i0 + j < NP && u[j] >= lo;
+    int total;
+    int at = n + block_scan(sm, m, phase, total);
+#pragma unroll
+    for (int j = 0; j < PER; ++j)
+      if (i0 + j < NP && u[j] >= lo) {
+        if (at < CAP) {
+          sm.key[at] = u[j];
+          sm.idx[at] = i0 + j;
+        }
+        at++;
+      }
+    n += total;
+    if (n > CAP) break;  // the whole row's then: no need to count on
+  }
+  __syncthreads();  // the candidates other threads stored
+  const bool fits = n >= K && n <= CAP;
+  if (fits) select(sm, Cand{&sm, n}, K, out, phase);
+  else select(sm, src, K, out, phase);
+  if (PATH != nullptr && threadIdx.x == 0) PATH[r] = fits ? 1 : 2;
+}
diff --git a/zig/kernels/cuda/glm_index_scores.cuh b/zig/kernels/cuda/glm_index_scores.cuh
new file mode 100644
index 0000000000..f6b83e7842
--- /dev/null
+++ b/zig/kernels/cuda/glm_index_scores.cuh
@@ -0,0 +1,130 @@
+// GLM-5.3-Flash DSA indexer, a prompt chunk's pool scores: tf_glm_scores' arithmetic with each key tile staged once.
+#pragma once
+
+#include <cuda_bf16.h>
+#include <math.h>
+#include <stdint.h>
+
+namespace glm_index {
+
+constexpr int D = 128, BP = 64, RB = 8, KQ = D / 16;  // key width, pools a tile, rows a CTA (a warp each), k16 steps
+
+#ifndef GLM_INDEX_HOST  // tools/glm_index_emulate supplies these on the host
+__device__ __forceinline__ void cp16(void* dst, const void* src, bool ok) {
+  const unsigned d = static_cast<unsigned>(__cvta_generic_to_shared(dst));
+  asm volatile("cp.async.cg.shared.global [%0], [%1], 16, %2;\n" ::"r"(d), "l"(src), "r"(ok ? 16 : 0));
+}
+__device__ __forceinline__ void cp_commit() { asm volatile("cp.async.commit_group;\n" ::); }
+__device__ __forceinline__ void cp_wait() { asm volatile("cp.async.wait_group 0;\n" ::); }
+
+__device__ __forceinline__ void ldm4(uint32_t* r, const void* p) {
+  const unsigned a = static_cast<unsigned>(__cvta_generic_to_shared(p));
+  asm volatile("ldmatrix.sync.aligned.m8n8.x4.shared.b16 {%0, %1, %2, %3}, [%4];\n"
+               : "=r"(r[0]), "=r"(r[1]), "=r"(r[2]), "=r"(r[3])
+               : "r"(a));
+}
+
+__device__ __forceinline__ void mma16816(float* c, const uint32_t* a, const uint32_t* b) {
+  asm volatile(
+      "mma.sync.aligned.m16n8k16.row.col.f32.bf16.bf16.f32 {%0, %1, %2, %3}, {%4, %5, %6, %7}, {%8, %9}, "
+      "{%0, %1, %2, %3};"
+      : "+f"(c[0]), "+f"(c[1]), "+f"(c[2]), "+f"(c[3])
+      : "r"(a[0]), "r"(a[1]), "r"(a[2]), "r"(a[3]), "r"(b[0]), "r"(b[1]));
+}
+#endif
+
+// Element offset of pool p's key: in order (pt null), or row p % 128 of page pt[p / 128] (glm_pages::at<128>).
+__device__ __forceinline__ int64_t key_at(const int* pt, int64_t page_elems, int64_t p) {
+  return pt ? (int64_t)pt[p / 128] * page_elems + (p % 128) * D : p * D;
+}
+
+// Logical 16-byte chunk u of key row `row` (0..63) in the tile: the 16-row groups' layout of tf_glm_scores.
+__device__ __forceinline__ int tile_at(int row, int u) { return row * D + ((u ^ (row & 7)) << 3); }
+
+}  // namespace glm_index
+
+// tf_glm_scores bit for bit, 8 rows (a warp each) sharing each 64-pool key tile; grid (cdiv(rows, 8), cdiv(nb, loop)).
+extern "C" __global__ void __launch_bounds__(glm_index::RB * 32)
+    tf_glm_index_scores(const __nv_bfloat16* qi, const __nv_bfloat16* wts, int64_t w_stride, const __nv_bfloat16* pk,
+                        float* out, const int* pos, int64_t rows, int64_t np, int64_t nb, float scale, float wscale,
+                        int64_t heads, int loop, const int* pt, int64_t ps) {
+  using namespace glm_index;
+  __shared__ __align__(16) __nv_bfloat16 tile[2][BP * D];  // two key tiles
+  const int lane = threadIdx.x & 31, warp = threadIdx.x >> 5, g = lane >> 2, q = lane & 3;
+  const int64_t r = (int64_t)blockIdx.x * RB + warp, P = pos[0];
+  const bool live = r < rows;
+  const int64_t npool = live ? (P + r + 1) / 4 : 0;
+  const int64_t last = min(rows, (int64_t)(blockIdx.x + 1) * RB) - 1, nmax = (P + last + 1) / 4;  // the CTA's largest
+  // this row's query fragments for the 8 k16 steps and 2 head tiles (ldmatrix's A layout), and a lane's 4 head weights
+  uint32_t aq[KQ][2][4];
+#pragma unroll
+  for (int kq = 0; kq < KQ; ++kq)
+#pragma unroll
+    for (int i = 0; i < 2; ++i) {
+      const int h0 = i * 16 + g, h1 = h0 + 8;
+      const __nv_bfloat16* b0 = qi + ((live ? r : 0) * heads + (h0 < heads ? h0 : 0)) * D + kq * 16 + 2 * q;
+      const __nv_bfloat16* b1 = qi + ((live ? r : 0) * heads + (h1 < heads ? h1 : 0)) * D + kq * 16 + 2 * q;
+      const bool o0 = live && h0 < heads, o1 = live && h1 < heads;
+      aq[kq][i][0] = o0 ? *reinterpret_cast<const uint32_t*>(b0) : 0u;
+      aq[kq][i][1] = o1 ? *reinterpret_cast<const uint32_t*>(b1) : 0u;
+      aq[kq][i][2] = o0 ? *reinterpret_cast<const uint32_t*>(b0 + 8) : 0u;
+      aq[kq][i][3] = o1 ? *reinterpret_cast<const uint32_t*>(b1 + 8) : 0u;
+    }
+  float wv[4];
+  for (int i = 0; i < 4; ++i) {
+    const int64_t hh = g + 8 * i;
+    wv[i] = live && hh < heads ? __fmul_rn(__bfloat162float(wts[r * w_stride + hh]), wscale) : 0.0f;
+  }
+  // The CTA's 256 threads copy tile pb's keys (pools past the CTA's visible ones read as 0).
+  auto stage = [&](int buf, int64_t pb) {
+    for (int idx = threadIdx.x; idx < BP * (D / 8); idx += RB * 32) {
+      const int row = idx / (D / 8), u = idx % (D / 8);
+      const int64_t p = pb * BP + row;
+      const bool ok = p < nmax;
+      cp16(&tile[buf][tile_at(row, u)], pk + key_at(pt, ps / 2, ok ? p : 0) + u * 8, ok);
+    }
+  };
+  const int64_t first = (int64_t)blockIdx.y * loop;
+  auto any = [&](int j) { return j < loop && first + j < nb && (first + j) * BP < nmax; };
+  if (any(0)) stage(0, first);
+  cp_commit();
+  for (int j = 0; j < loop; ++j) {
+    const int64_t pb = first + j, p0 = pb * BP;
+    if (pb >= nb) break;
+    cp_wait();
+    __syncthreads();  // tile j is in, and every warp is done with tile j - 1
+    if (any(j + 1)) stage((j + 1) & 1, pb + 1);
+    cp_commit();
+    if (!live) continue;
+    if (p0 >= npool) {  // the row sees none of these pools: -inf
+      for (int e = lane; e < BP; e += 32)
+        if (p0 + e < np) out[r * np + p0 + e] = -INFINITY;
+      continue;
+    }
+    for (int gi = 0; gi < BP / 16; ++gi) {
+      float acc[2][2][4] = {};
+      const __nv_bfloat16* keys = tile[j & 1] + gi * 16 * D;
+#pragma unroll
+      for (int kq = 0; kq < KQ; ++kq) {
+        uint32_t b[4];
+        ldm4(b, keys + tile_at((lane >> 4) * 8 + (lane & 7), kq * 2 + ((lane >> 3) & 1)));
+        for (int nt = 0; nt < 2; ++nt)
+          for (int i = 0; i < 2; ++i) mma16816(acc[i][nt], aq[kq][i], b + 2 * nt);
+      }
+      for (int nt = 0; nt < 2; ++nt)
+        for (int e = 0; e < 2; ++e) {
+          const float v0 = fmaxf(__fmul_rn(acc[0][nt][e], scale), 0.0f),
+                      v1 = fmaxf(__fmul_rn(acc[0][nt][e + 2], scale), 0.0f);
+          const float v2 = fmaxf(__fmul_rn(acc[1][nt][e], scale), 0.0f),
+                      v3 = fmaxf(__fmul_rn(acc[1][nt][e + 2], scale), 0.0f);
+          float s = __fmul_rn(wv[1], v1);
+          s = __fmaf_rn(wv[0], v0, s);
+          s = __fmaf_rn(wv[2], v2, s);
+          s = __fmaf_rn(wv[3], v3, s);
+          for (int o = 16; o >= 4; o >>= 1) s = s + __shfl_xor_sync(0xffffffffu, s, o);
+          const int64_t p = p0 + gi * 16 + nt * 8 + q * 2 + e;
+          if (g == 0 && p < np) out[r * np + p] = p < npool ? s : -INFINITY;
+        }
+    }
+  }
+}
diff --git a/zig/src/core/contract/glm_cuda_records.zig b/zig/src/core/contract/glm_cuda_records.zig
index 1f0dee987f..540ca63d80 100644
--- a/zig/src/core/contract/glm_cuda_records.zig
+++ b/zig/src/core/contract/glm_cuda_records.zig
@@ -41,6 +41,8 @@ pub const cuda = [_]registry.CudaRecord{
     .{ .launch = .{ .name = "glm_glue tf_glm_lat_write", .symbol = "tf_glm_lat_write", .image = "glm_glue", .threads = 128, .min_sm = 80 }, .usage = .{ .{ .regs = 13, .static = 0, .bound = 0 }, .{ .regs = 13, .static = 0, .bound = 0 }, .{ .regs = 13, .static = 0, .bound = 0 }, .{ .regs = 36, .static = 0, .bound = 0 }, .{ .regs = 36, .static = 0, .bound = 0 } } },
     .{ .launch = .{ .name = "glm_glue tf_glm_pool_keys", .symbol = "tf_glm_pool_keys", .image = "glm_glue", .threads = 32, .min_sm = 80 }, .usage = .{ .{ .regs = 25, .static = 0, .bound = 32 }, .{ .regs = 26, .static = 0, .bound = 32 }, .{ .regs = 26, .static = 0, .bound = 32 }, .{ .regs = 24, .static = 0, .bound = 32 }, .{ .regs = 24, .static = 0, .bound = 32 } } },
     .{ .launch = .{ .name = "glm_glue tf_glm_scores", .symbol = "tf_glm_scores", .image = "glm_glue", .threads = 128, .min_sm = 80 }, .usage = .{ .{ .regs = 126, .static = 32768, .bound = 128 }, .{ .regs = 131, .static = 32768, .bound = 128 }, .{ .regs = 131, .static = 32768, .bound = 128 }, .{ .regs = 155, .static = 32768, .bound = 128 }, .{ .regs = 155, .static = 32768, .bound = 128 } } },
+    .{ .launch = .{ .name = "glm_index tf_glm_index_scores", .symbol = "tf_glm_index_scores", .image = "glm_index", .threads = 256, .min_sm = 80 }, .usage = .{ .{ .regs = 128, .static = 32768, .bound = 256 }, .{ .regs = 128, .static = 32768, .bound = 256 }, .{ .regs = 128, .static = 32768, .bound = 256 }, .{ .regs = 124, .static = 32768, .bound = 256 }, .{ .regs = 124, .static = 32768, .bound = 256 } } },
+    .{ .launch = .{ .name = "glm_index tf_glm_index_pools", .symbol = "tf_glm_index_pools", .image = "glm_index", .threads = 256, .min_sm = 80 }, .usage = .{ .{ .regs = 62, .static = 33864, .bound = 256 }, .{ .regs = 64, .static = 33864, .bound = 256 }, .{ .regs = 64, .static = 33864, .bound = 256 }, .{ .regs = 59, .static = 33864, .bound = 256 }, .{ .regs = 59, .static = 33864, .bound = 256 } } },
     .{ .launch = .{ .name = "glm_glue tf_glm_chunks_16", .symbol = "tf_glm_chunks_16", .image = "glm_glue", .threads = 256, .dynamic = 83328, .opt_in = true, .min_sm = 80 }, .usage = .{ .{ .regs = 64, .static = 0, .bound = 256 }, .{ .regs = 64, .static = 0, .bound = 256 }, .{ .regs = 64, .static = 0, .bound = 256 }, .{ .regs = 123, .static = 0, .bound = 256 }, .{ .regs = 123, .static = 0, .bound = 256 } } },
     .{ .launch = .{ .name = "glm_glue tf_glm_chunks_32", .symbol = "tf_glm_chunks_32", .image = "glm_glue", .threads = 256, .dynamic = 101120, .opt_in = true, .min_sm = 80 }, .usage = .{ .{ .regs = 150, .static = 0, .bound = 256 }, .{ .regs = 150, .static = 0, .bound = 256 }, .{ .regs = 150, .static = 0, .bound = 256 }, .{ .regs = 130, .static = 0, .bound = 256 }, .{ .regs = 130, .static = 0, .bound = 256 } } },
     .{ .launch = .{ .name = "glm_ops tf_glm_dsa_tokens", .symbol = "tf_glm_dsa_tokens", .image = "glm_ops", .threads = 256, .min_sm = 80 }, .usage = .{ .{ .regs = 25, .static = 0, .bound = 0 }, .{ .regs = 25, .static = 0, .bound = 0 }, .{ .regs = 25, .static = 0, .bound = 0 }, .{ .regs = 24, .static = 0, .bound = 0 }, .{ .regs = 24, .static = 0, .bound = 0 } } },
@@ -60,6 +62,8 @@ pub const cuda_132 = [_]registry.Override{
 pub const cuda_133 = [_]registry.Override{
     .{ .name = "glm_glue tf_glm_lat_write", .usage = .{ .{ .regs = 13, .static = 0, .bound = 0 }, .{ .regs = 13, .static = 0, .bound = 0 }, .{ .regs = 13, .static = 0, .bound = 0 }, .{ .regs = 34, .static = 0, .bound = 0 }, .{ .regs = 34, .static = 0, .bound = 0 } } },
     .{ .name = "glm_glue tf_glm_scores", .usage = .{ .{ .regs = 122, .static = 32768, .bound = 128 }, .{ .regs = 129, .static = 32768, .bound = 128 }, .{ .regs = 129, .static = 32768, .bound = 128 }, .{ .regs = 157, .static = 32768, .bound = 128 }, .{ .regs = 157, .static = 32768, .bound = 128 } } },
+    .{ .name = "glm_index tf_glm_index_scores", .usage = .{ .{ .regs = 128, .static = 32768, .bound = 256 }, .{ .regs = 128, .static = 32768, .bound = 256 }, .{ .regs = 128, .static = 32768, .bound = 256 }, .{ .regs = 126, .static = 32768, .bound = 256 }, .{ .regs = 126, .static = 32768, .bound = 256 } } },
+    .{ .name = "glm_index tf_glm_index_pools", .usage = .{ .{ .regs = 62, .static = 33864, .bound = 256 }, .{ .regs = 64, .static = 33864, .bound = 256 }, .{ .regs = 64, .static = 33864, .bound = 256 }, .{ .regs = 54, .static = 33864, .bound = 256 }, .{ .regs = 54, .static = 33864, .bound = 256 } } },
     .{ .name = "glm_glue tf_glm_chunks_32", .usage = .{ .{ .regs = 150, .static = 0, .bound = 256 }, .{ .regs = 150, .static = 0, .bound = 256 }, .{ .regs = 150, .static = 0, .bound = 256 }, .{ .regs = 128, .static = 0, .bound = 256 }, .{ .regs = 128, .static = 0, .bound = 256 } } },
     .{ .name = "glm_ops tf_glm_dsa_tokens", .usage = .{ .{ .regs = 25, .static = 0, .bound = 0 }, .{ .regs = 25, .static = 0, .bound = 0 }, .{ .regs = 25, .static = 0, .bound = 0 }, .{ .regs = 22, .static = 0, .bound = 0 }, .{ .regs = 22, .static = 0, .bound = 0 } } },
 };
diff --git a/zig/src/cuda/kernels.zig b/zig/src/cuda/kernels.zig
index 6037e7bb08..33a3bbfe55 100644
--- a/zig/src/cuda/kernels.zig
+++ b/zig/src/cuda/kernels.zig
@@ -68,6 +68,7 @@ pub const dense_rows: []const u8 = if (available) &Blob("fatbin_dense_rows").byt
 pub const modelopt_pack: []const u8 = if (available) &Blob("fatbin_modelopt_pack").bytes else &.{};
 pub const glm_glue: []const u8 = if (available) &Blob("fatbin_glm_glue").bytes else &.{};
 pub const glm_ops: []const u8 = if (available) &Blob("fatbin_glm_ops").bytes else &.{};
+pub const glm_index: []const u8 = if (available) &Blob("fatbin_glm_index").bytes else &.{};
 pub const glm_quant: []const u8 = if (available) &Blob("fatbin_glm_quant").bytes else &.{};
 pub const glm_vision: []const u8 = if (available) &Blob("fatbin_glm_vision").bytes else &.{};
 pub const glm_mm: []const u8 = if (available) &Blob("fatbin_glm_mm").bytes else &.{};
diff --git a/zig/src/families/glm/cuda.zig b/zig/src/families/glm/cuda.zig
index 5557aa292b..6a11904a3c 100644
--- a/zig/src/families/glm/cuda.zig
+++ b/zig/src/families/glm/cuda.zig
@@ -13,6 +13,7 @@ pub const kda = @import("cuda_kda.zig");
 pub const b16 = @import("cuda_b16.zig");
 pub const latent = @import("cuda_latent.zig");
 pub const sparse = @import("cuda_sparse.zig");
+pub const index = @import("cuda_index.zig");
 pub const moe = @import("cuda_moe.zig");
 pub const kda_block = @import("cuda_kda_block.zig");
 pub const dsa_block = @import("cuda_dsa_block.zig");
@@ -56,6 +57,7 @@ test {
     _ = b16;
     _ = latent;
     _ = sparse;
+    _ = index;
     _ = moe;
     _ = kda_block;
     _ = dsa_block;
diff --git a/zig/src/families/glm/cuda_forward.zig b/zig/src/families/glm/cuda_forward.zig
index 7381abe804..405097e687 100644
--- a/zig/src/families/glm/cuda_forward.zig
+++ b/zig/src/families/glm/cuda_forward.zig
@@ -24,6 +24,7 @@ const prof_mod = @import("cuda_prof.zig");
 const images_mod = @import("cuda_images.zig");
 const exl3 = @import("cuda_exl3.zig");
 const shared = @import("cuda_shared.zig");
+const index_mod = @import("cuda_index.zig");
 
 pub const prefill_rows = state.prefill_rows;
 pub const window_rows = state.window_rows;
@@ -40,6 +41,7 @@ pub const Modules = struct {
     experts: experts.Kernels,
     exl3: exl3.Kernels, // an EXL3 checkpoint's routed experts
     mm: mm_mod.Kernels, // glm_mm: the bf16 linears on tiled weights
+    index: ?index_mod.Kernels, // TF_GLM_INDEX_CUDA: a prompt chunk's indexer scores and pools (glm_index)
     tiled: bool = false, // eligible bf16 linears are tiled; KDA inputs with uneven strides retain their row-major path
     packs: ?*const mm_mod.Packed = null, // the linears packed to 4, 6 or 8 bits (TF_GLM_DENSE_BITS, _HEAD_BITS)
     copy_rows: cuda.Function,
@@ -73,6 +75,13 @@ pub const Modules = struct {
         m.exl3 = try exl3.Kernels.load(d);
         errdefer m.exl3.unload();
         m.mm = try mm_mod.Kernels.load(m.mods[7]);
+        m.index = null;
+        errdefer if (m.index) |*x| x.unload();
+        const mode = index_mod.Mode.read() catch return error.BadIndexSetting;
+        if (mode != .off) m.index = index_mod.Kernels.load(d, k.glm_index, mode) catch |err| blk: {
+            std.log.warn("TF_GLM_INDEX_CUDA: the indexer kernels are not used ({t}); prompt chunks use tf_glm_scores", .{err});
+            break :blk null;
+        };
         m.tiled = false;
         m.packs = null; // `m` starts undefined: the field default does not apply
         return m;
@@ -88,6 +97,7 @@ pub const Modules = struct {
     }
 
     pub fn deinit(m: *Modules) void {
+        if (m.index) |*x| x.unload();
         m.exl3.unload();
         for (&m.mods) |*x| x.unload();
         m.lane.unload();
@@ -184,7 +194,7 @@ pub const Forward = struct {
     taps: ?taps_mod.Taps = null, // DFlash2 (--drafter): the tapped layers' rows, after each such layer
 
     pub fn set(f: *const Forward) glue.Set {
-        return .{ .glue = &f.m.glue, .mm = if (f.m.tiled) &f.m.mm else null, .packs = f.m.packs, .prompt = f.ranks[0].b.prefill, .linear = .{ .lane = &f.m.lane }, .copy_rows = f.m.copy_rows };
+        return .{ .glue = &f.m.glue, .index = if (f.m.index) |*x| x else null, .mm = if (f.m.tiled) &f.m.mm else null, .packs = f.m.packs, .prompt = f.ranks[0].b.prefill, .linear = .{ .lane = &f.m.lane }, .copy_rows = f.m.copy_rows };
     }
 
     /// reduce.sums' gather: every rank's partial of `rows` rows into every rank's gathered buffer, rank order.
diff --git a/zig/src/families/glm/cuda_glue.zig b/zig/src/families/glm/cuda_glue.zig
index 4d0343c348..5f044f0b20 100644
--- a/zig/src/families/glm/cuda_glue.zig
+++ b/zig/src/families/glm/cuda_glue.zig
@@ -3,9 +3,10 @@
 const std = @import("std");
 const cuda = @import("cuda");
 const mm_mod = @import("cuda_mm.zig");
+const index_mod = @import("cuda_index.zig");
 
 /// What the forward's kernels launch from: our glue's, glm_mm's when tiled, qlinear's lane for FP8 kinds, copy_rows.
-pub const Set = struct { glue: *const Kernels, mm: ?*const mm_mod.Kernels = null, packs: ?*const mm_mod.Packed = null, prompt: bool = false, linear: cuda.qlinear.Linear = .{}, copy_rows: ?cuda.Function = null };
+pub const Set = struct { glue: *const Kernels, index: ?*const index_mod.Kernels = null, mm: ?*const mm_mod.Kernels = null, packs: ?*const mm_mod.Packed = null, prompt: bool = false, linear: cuda.qlinear.Linear = .{}, copy_rows: ?cuda.Function = null };
 
 /// A hidden-size norm's weight: bf16, or float32 (FP4-Local's) read as float32 by the _f32w kernels.
 pub const Norm = union(enum) { bf16: u64, f32: u64 };
diff --git a/zig/src/families/glm/cuda_index.zig b/zig/src/families/glm/cuda_index.zig
new file mode 100644
index 0000000000..5c4ad55bb2
--- /dev/null
+++ b/zig/src/families/glm/cuda_index.zig
@@ -0,0 +1,131 @@
+//! A prompt chunk's DSA indexer in CUDA (glm_index.cu): the bytes of tf_glm_scores and tf_glm_select_rows.
+
+const std = @import("std");
+const cuda = @import("cuda");
+
+const sparse = @import("cuda_sparse.zig");
+
+pub const topk = sparse.topk_pools;
+/// glm_index_scores: rows a CTA (a warp each), 64-pool tiles a CTA, threads a CTA.
+pub const score_rows = 8;
+pub const score_tiles = 16;
+pub const threads = score_rows * 32;
+/// glm_index_pools' block: glm_index::T.
+pub const pool_threads = 256;
+/// The rows a launch must have to take it: where tf_glm_scores' `loop` tiling starts.
+pub const rows_from = sparse.score_loop_from;
+pub const cap_max = 4096;
+
+/// TF_GLM_INDEX_CUDA: unset or 0 (today's launches), 1 (both kernels), scores (only the scores, to bisect).
+pub const Mode = enum {
+    off,
+    both,
+    scores,
+
+    pub fn parse(text: ?[]const u8) !Mode {
+        const t = text orelse return .off;
+        if (t.len == 0 or std.mem.eql(u8, t, "0")) return .off;
+        if (std.mem.eql(u8, t, "1")) return .both;
+        if (std.mem.eql(u8, t, "scores")) return .scores;
+        return error.BadIndexMode;
+    }
+
+    pub fn read() !Mode {
+        const v = std.c.getenv("TF_GLM_INDEX_CUDA") orelse return .off;
+        return parse(std.mem.span(v));
+    }
+};
+
+/// glm_index_scores' grid: rows in CTAs of 8, the `nb` pool blocks in runs of score_tiles.
+pub fn scoresGrid(rows: usize, nb: usize) [2]u32 {
+    return .{ @intCast(cdiv(rows, score_rows)), @intCast(@max(1, cdiv(nb, score_tiles))) };
+}
+
+/// The candidates a row of the pools kernel keeps and the floor's target (speed only; the pools are the same).
+pub const Tile = struct { cap: usize, target: usize };
+
+pub fn tile(np: usize) Tile {
+    if (np <= 12288) return .{ .cap = 1024, .target = 724 };
+    if (np <= 32768) return .{ .cap = 2048, .target = 1024 };
+    return .{ .cap = cap_max, .target = 2048 };
+}
+
+pub const Kernels = struct {
+    module: cuda.Module,
+    scores_fn: cuda.Function,
+    pools_fn: cuda.Function,
+    /// Whether the pools kernel replaces tf_glm_select_rows (Mode.both).
+    select: bool,
+
+    /// `image`: cuda.kernels.glm_index.
+    pub fn load(d: *const cuda.Driver, image: []const u8, mode: Mode) !Kernels {
+        if (image.len == 0) return error.NoIndexKernels;
+        var m = try cuda.Module.load(d, image);
+        errdefer m.unload();
+        return .{ .module = m, .scores_fn = try m.function("tf_glm_index_scores"), .pools_fn = try m.function("tf_glm_index_pools"), .select = mode == .both };
+    }
+
+    pub fn unload(k: *Kernels) void {
+        k.module.unload();
+    }
+
+    /// Whether `rows` rows over `np` pools with this key cache layout go through the kernels.
+    pub fn takes(k: *const Kernels, rows: usize, np: usize, pk: u64, ps: u64) bool {
+        _ = k;
+        return rows >= rows_from and np >= topk and pk % 16 == 0 and ps % 16 == 0;
+    }
+
+    /// tf_glm_scores' launch for `rows` rows at `at` into out [rows][np]; the first `nb` 64-pool blocks are written.
+    pub fn scores(k: *const Kernels, s: cuda.Stream, qi: u64, wts: u64, w_stride: usize, pk: u64, out: u64, at: u64, rows: usize, np: usize, nb: usize, scale: f32, wscale: f32, heads: usize, pt: u64, ps: u64) !void {
+        const g = scoresGrid(rows, nb);
+        var a: cuda.Args = .{};
+        inline for (.{ qi, wts, @as(i64, @intCast(w_stride)), pk, out, at, @as(i64, @intCast(rows)), @as(i64, @intCast(np)), @as(i64, @intCast(nb)), scale, wscale, @as(i64, @intCast(heads)), @as(c_int, score_tiles), pt, @as(i64, @intCast(ps)) }) |v| a.add(v);
+        try cuda.launch.launch(k.scores_fn, .{ .grid = .{ .x = g[0], .y = g[1] }, .block = .{ .x = threads } }, s, &a);
+    }
+
+    /// tf_glm_select_rows' launch: each row's 512 best pools of scores [rows][np] ascending into out [rows][512] int64.
+    pub fn pools(k: *const Kernels, s: cuda.Stream, scores_ptr: u64, out: u64, at: u64, rows: usize, np: usize, path: u64) !void {
+        const t = tile(np);
+        var a: cuda.Args = .{};
+        inline for (.{ scores_ptr, out, @as(i64, @intCast(np)), at, @as(c_int, topk), @as(c_int, @intCast(t.cap)), @as(c_int, @intCast(t.target)), path }) |v| a.add(v);
+        try cuda.launch.launch(k.pools_fn, .{ .grid = .{ .x = @intCast(rows) }, .block = .{ .x = pool_threads } }, s, &a);
+    }
+};
+
+fn cdiv(a: usize, b: usize) usize {
+    return (a + b - 1) / b;
+}
+
+test "GLM index: TF_GLM_INDEX_CUDA is off unless 1 or scores, and a typo is refused" {
+    try std.testing.expectEqual(Mode.off, try Mode.parse(null));
+    try std.testing.expectEqual(Mode.off, try Mode.parse(""));
+    try std.testing.expectEqual(Mode.off, try Mode.parse("0"));
+    try std.testing.expectEqual(Mode.both, try Mode.parse("1"));
+    try std.testing.expectEqual(Mode.scores, try Mode.parse("scores"));
+    try std.testing.expectError(error.BadIndexMode, Mode.parse("on"));
+}
+
+test "GLM index: the scores grid covers 8 rows a CTA and 16 pool blocks a CTA" {
+    // a 512-row chunk over 147,456 pools: 2,304 blocks
+    try std.testing.expectEqual([2]u32{ 64, 144 }, scoresGrid(512, 2304));
+    try std.testing.expectEqual([2]u32{ 9, 1 }, scoresGrid(65, 16));
+    try std.testing.expectEqual([2]u32{ 8, 2 }, scoresGrid(64, 17));
+    try std.testing.expectEqual([2]u32{ 1, 1 }, scoresGrid(1, 0));
+}
+
+test "GLM index: the candidate cap follows the pool count and stays in shared memory" {
+    try std.testing.expectEqual(Tile{ .cap = 1024, .target = 724 }, tile(12288));
+    try std.testing.expectEqual(Tile{ .cap = 2048, .target = 1024 }, tile(12289));
+    try std.testing.expectEqual(Tile{ .cap = 2048, .target = 1024 }, tile(32768));
+    try std.testing.expectEqual(Tile{ .cap = 4096, .target = 2048 }, tile(32769));
+    try std.testing.expect(tile(1 << 30).cap <= cap_max);
+}
+
+test "GLM index: only prompt-sized launches on aligned key caches take the kernels" {
+    const k: Kernels = undefined;
+    try std.testing.expect(k.takes(rows_from, 1024, 0x1000, 32768));
+    try std.testing.expect(!k.takes(rows_from - 1, 1024, 0x1000, 32768));
+    try std.testing.expect(!k.takes(512, 1024, 0x1008, 32768));
+    try std.testing.expect(!k.takes(512, 1024, 0x1000, 32760));
+    try std.testing.expect(!k.takes(512, topk - 1, 0x1000, 32768));
+}
diff --git a/zig/src/families/glm/cuda_prompt_logprobs.zig b/zig/src/families/glm/cuda_prompt_logprobs.zig
index c1bcecd738..3694a14358 100644
--- a/zig/src/families/glm/cuda_prompt_logprobs.zig
+++ b/zig/src/families/glm/cuda_prompt_logprobs.zig
@@ -113,7 +113,7 @@ test "the head's columns at TP=2 and TP=3 take no K slices: a pass of any rows n
     };
 }
 
-/// The projection policy used by the prompt-row head passes; an MXFP8 head (FP4-Local's) runs on the lane.
+/// The head passes' projection policy (no indexer launches, so no .index); an MXFP8 head runs on the lane.
 pub fn projectionSet(e: anytype) @import("cuda_glue.zig").Set {
     return .{ .glue = &e.m.glue, .mm = if (e.m.tiled) &e.m.mm else null, .packs = e.m.packs, .prompt = true, .linear = .{ .lane = &e.m.lane } };
 }
diff --git a/zig/src/families/glm/cuda_sparse.zig b/zig/src/families/glm/cuda_sparse.zig
index c09bf39eab..cccd010e59 100644
--- a/zig/src/families/glm/cuda_sparse.zig
+++ b/zig/src/families/glm/cuda_sparse.zig
@@ -73,13 +73,20 @@ pub fn selectTokensProfile(set: glue.Set, ops: cuda.Module, stream: cuda.Stream,
             at = s.at;
         }
         const nb = cdiv(scoreColumns(host_pos + a, n, np_max), 64);
+        const scale: f32 = @floatCast(std.math.pow(f64, @floatFromInt(d), -0.5));
         var sa: cuda.Args = .{};
-        inline for (.{ qi + a * heads * d * 2, wts + a * w_stride * 2, @as(i64, @intCast(w_stride)), pk, s.scores, at, @as(i64, @intCast(np_max)), @as(i64, @intCast(nb)), @as(f32, @floatCast(std.math.pow(f64, @floatFromInt(d), -0.5))), headScale(heads), @as(i64, @intCast(heads)), @as(c_int, @intCast(loop)), t.pt, @as(i64, @intCast(t.ps)) }) |v| sa.add(v);
+        inline for (.{ qi + a * heads * d * 2, wts + a * w_stride * 2, @as(i64, @intCast(w_stride)), pk, s.scores, at, @as(i64, @intCast(np_max)), @as(i64, @intCast(nb)), scale, headScale(heads), @as(i64, @intCast(heads)), @as(c_int, @intCast(loop)), t.pt, @as(i64, @intCast(t.ps)) }) |v| sa.add(v);
         if (d != 128 or heads > 32) return error.IndexShape;
+        const quick = if (set.index) |k| if (k.takes(rows, np_max, pk, t.ps)) k else null else null;
         var stamp = prof.span(profile, stream, .index_score);
-        try cuda.launch.launch(set.glue.scores, .{ .grid = .{ .x = @intCast(n), .y = @intCast(cdiv(nb, loop)) }, .block = .{ .x = 128 } }, stream, &sa);
+        if (quick) |k| try k.scores(stream, qi + a * heads * d * 2, wts + a * w_stride * 2, w_stride, pk, s.scores, at, n, np_max, nb, scale, headScale(heads), heads, t.pt, t.ps) else try cuda.launch.launch(set.glue.scores, .{ .grid = .{ .x = @intCast(n), .y = @intCast(cdiv(nb, loop)) }, .block = .{ .x = 128 } }, stream, &sa);
         stamp.end();
         stamp = prof.span(profile, stream, .index_select);
+        if (quick) |k| if (k.select) {
+            try k.pools(stream, s.scores, s.pools + a * topk_pools * 8, at, n, np_max, 0);
+            stamp.end();
+            continue;
+        };
         var sel: cuda.Args = .{};
         inline for (.{ s.scores, s.pools + a * topk_pools * 8, @as(i64, @intCast(np_max)), at, @as(c_int, topk_pools) }) |v| sel.add(v);
         try cuda.launch.launch(set.glue.select_rows, .{ .grid = .{ .x = @intCast(n) }, .block = .{ .x = 128 } }, stream, &sel);
diff --git a/zig/tests/cuda/contract_glm_index.zig b/zig/tests/cuda/contract_glm_index.zig
index 12ded25c8f..3c5736ac1e 100644
--- a/zig/tests/cuda/contract_glm_index.zig
+++ b/zig/tests/cuda/contract_glm_index.zig
@@ -41,6 +41,9 @@ const sumIn = pages.sumIn;
 pub const Scores = struct {
     rig: rig_.Rig,
     f: cuda.Function,
+    /// tf_glm_index_scores, which must equal `f`'s bytes (`f` is then the baseline it replaces).
+    fast: ?cuda.Function = null,
+    base: u64 = 0,
     wrong_table: bool = false,
     shape: c.Shape = undefined,
     q: []u16 = &.{},
@@ -73,6 +76,10 @@ pub const Scores = struct {
 
     pub fn kernel(self: *Scores) g.Kernel {
         var declared = spec(.scores);
+        if (self.fast != null) {
+            declared.name = "cuda glm tf_glm_index_scores on pages out of order";
+            declared.replaces = "tf_glm_scores";
+        }
         if (self.wrong_table) declared.name = "control: GLM paged scores use a swapped table";
         return .{ .spec = declared, .ctx = self, .vt = &vt, .control = if (self.wrong_table) .shapes else null };
     }
@@ -125,17 +132,34 @@ pub const Scores = struct {
         self.pos = try self.rig.input(std.mem.asBytes(&p));
         self.pt = try pages.uploadTable(&self.rig, self.ids, self.wrong_table);
         self.out = try self.rig.output(shape.rows * shape.n * 4);
+        self.base = try self.rig.scratch(shape.rows * shape.n * 4);
         return true;
     }
 
     pub fn run(self: *Scores) !void {
+        try self.launch(self.fast orelse self.f, self.fast != null, self.out);
+        try self.rig.s.synchronize();
+    }
+
+    fn launch(self: *Scores, f: cuda.Function, fast: bool, out: u64) !void {
         const rows = self.shape.rows;
         const nb = self.shape.n / 64;
-        const loop: usize = if (rows >= sparse.score_loop_from) sparse.score_loop else 1;
+        const loop: usize = if (fast) glm.index.score_tiles else if (rows >= sparse.score_loop_from) sparse.score_loop else 1;
         var a: cuda.Args = .{};
-        inline for (.{ self.in[0], self.in[1], @as(i64, w_stride), self.in[2], self.out, self.pos, @as(i64, @intCast(self.shape.n)), @as(i64, @intCast(nb)), self.scale, self.wscale, @as(i64, heads), @as(c_int, @intCast(loop)), self.pt, @as(i64, pages.page_bytes) }) |v| a.add(v);
-        try cuda.launch.launch(self.f, .{ .grid = .{ .x = @intCast(rows), .y = @intCast((nb + loop - 1) / loop) }, .block = .{ .x = 128 } }, self.rig.s, &a);
-        try self.rig.s.synchronize();
+        if (fast) {
+            inline for (.{ self.in[0], self.in[1], @as(i64, w_stride), self.in[2], out, self.pos, @as(i64, @intCast(rows)), @as(i64, @intCast(self.shape.n)), @as(i64, @intCast(nb)), self.scale, self.wscale, @as(i64, heads), @as(c_int, @intCast(loop)), self.pt, @as(i64, pages.page_bytes) }) |v| a.add(v);
+            const g2 = glm.index.scoresGrid(rows, nb);
+            try cuda.launch.launch(f, .{ .grid = .{ .x = g2[0], .y = g2[1] }, .block = .{ .x = glm.index.threads } }, self.rig.s, &a);
+        } else {
+            inline for (.{ self.in[0], self.in[1], @as(i64, w_stride), self.in[2], out, self.pos, @as(i64, @intCast(self.shape.n)), @as(i64, @intCast(nb)), self.scale, self.wscale, @as(i64, heads), @as(c_int, @intCast(loop)), self.pt, @as(i64, pages.page_bytes) }) |v| a.add(v);
+            try cuda.launch.launch(f, .{ .grid = .{ .x = @intCast(rows), .y = @intCast((nb + loop - 1) / loop) }, .block = .{ .x = 128 } }, self.rig.s, &a);
+        }
+    }
+
+    /// tf_glm_scores on the same inputs, as cuda_sparse launches it.
+    pub fn replaced(self: *Scores) ![]const u8 {
+        try self.launch(self.f, false, self.base);
+        return self.rig.bytes(self.base, self.shape.rows * self.shape.n * 4);
     }
 
     pub fn count(self: *Scores) usize {
diff --git a/zig/tests/cuda/contract_glm_pools.zig b/zig/tests/cuda/contract_glm_pools.zig
new file mode 100644
index 0000000000..016e2cc8db
--- /dev/null
+++ b/zig/tests/cuda/contract_glm_pools.zig
@@ -0,0 +1,185 @@
+//! GLM's prompt-chunk select under the contract on CUDA: tf_glm_index_pools, bit-equal to tf_glm_select_rows.
+const std = @import("std");
+const cuda = @import("cuda");
+const core = @import("core");
+const glm = @import("glm_cuda");
+const rig_ = @import("contract_rig.zig");
+
+const contract = core.contract;
+const c = contract.spec;
+const g = contract.gates;
+const sparse = glm.sparse;
+
+fn noSums(_: c.Shape) usize {
+    return 0;
+}
+
+const rule = "tf_glm_select_rows' bytes, and the host's picks: order keys descending, -0 as +0 and any NaN above +inf " ++
+    "(the GPU's s + 0 gives its one positive NaN), ties to the lower pool, the K picks in pool order";
+
+/// A prompt block's rows over buckets that take the whole-row path (1k), the candidates (4k to 64k) and the fallback.
+pub fn spec() c.Spec {
+    const k = sparse.topk_pools;
+    const served = &[_]c.Shape{ .{ .rows = 64, .n = 4096, .k = k }, .{ .rows = 16, .n = 65536, .k = k } };
+    return .{ .name = "cuda glm select rows tf_glm_index_pools", .accum = .none, .depth = noSums, .out = .raw, .replaces = "tf_glm_select_rows", .nan = .{ .stated = rule }, .tiles = .{ .{}, .{ .size = 8192, .min = k }, .{} }, .served = served };
+}
+
+/// sparse._order_key as the GPU computes it: s + 0 turns -0 into +0 and every NaN into one positive NaN.
+fn orderKey(s: f32) u32 {
+    const b: i32 = if (std.math.isNan(s)) std.math.maxInt(i32) else @bitCast(s + 0.0);
+    return @bitCast(b ^ ((b >> 31) | std.math.minInt(i32)));
+}
+
+fn before(keys: []const u32, x: u32, y: u32) bool {
+    return keys[x] > keys[y] or (keys[x] == keys[y] and x < y);
+}
+
+/// Each row's K pools of an n-pool bucket (shape: rows x n x K); tf_glm_select_rows' bytes are its baseline.
+pub const Pools = struct {
+    rig: rig_.Rig,
+    fast: cuda.Function,
+    narrow: cuda.Function,
+    shape: c.Shape = undefined,
+    host: []f32 = &.{},
+    at: usize = 0,
+    in: u64 = 0,
+    pos: u64 = 0,
+    out: u64 = 0,
+    base: u64 = 0,
+    msg: [160]u8 = undefined,
+
+    pub fn deinit(self: *Pools) void {
+        self.rig.gpu.gpa.free(self.host);
+        self.rig.deinit();
+    }
+
+    pub fn kernel(self: *Pools) g.Kernel {
+        return .{ .spec = spec(), .ctx = self, .vt = &vt };
+    }
+
+    const vt = g.hooks(Pools);
+
+    fn visible(self: *const Pools, r: usize) usize {
+        return (self.at + r + 1) / sparse.pool;
+    }
+
+    fn outBytes(self: *const Pools) usize {
+        return self.shape.rows * self.shape.k * 8;
+    }
+
+    /// Quarter steps in [-512, 512) with -0 among the zeros, -inf past the visible pools; non-finite: rows 0-3 plants.
+    pub fn prepare(self: *Pools, shape: c.Shape, fill: g.Fill, seed: u64) !bool {
+        if (shape.k != sparse.topk_pools or shape.n < shape.k or shape.rows == 0) return false;
+        self.rig.reset();
+        const a = self.rig.gpu.gpa;
+        a.free(self.host);
+        self.host = &.{};
+        self.shape = shape;
+        // a bucket's three quarters visible; past an edge's n every pool is (a limit of n), at 512 fewer than K are
+        self.at = if (std.math.isPowerOfTwo(shape.n)) 3 * shape.n else 4 * shape.n;
+        self.host = try a.alloc(f32, shape.rows * shape.n);
+        var prng = std.Random.DefaultPrng.init(seed);
+        const rng = prng.random();
+        for (0..shape.rows) |r| {
+            const vis = self.visible(r);
+            for (self.host[r * shape.n ..][0..shape.n], 0..) |*v, p| {
+                const level: f32 = @floatFromInt(@as(i32, @intCast(rng.uintLessThan(u32, 4096))) - 2048);
+                v.* = if (p >= vis) -std.math.inf(f32) else if (level == 0 and rng.boolean()) -0.0 else level * 0.25;
+            }
+        }
+        if (fill == .nonfinite) {
+            const neg_nan: f32 = @bitCast(@as(u32, 0xffc00000));
+            const put = [_]f32{ std.math.nan(f32), neg_nan, std.math.inf(f32), -std.math.inf(f32) };
+            for (put, 0..) |v, r| if (r < shape.rows) {
+                const step = @min(self.visible(r), shape.n) / 8;
+                for (0..8) |e| self.host[r * shape.n + e * step + r] = v;
+            };
+        }
+        const p: i32 = @intCast(self.at);
+        self.in = try self.rig.input(std.mem.sliceAsBytes(self.host));
+        self.pos = try self.rig.input(std.mem.asBytes(&p));
+        self.out = try self.rig.output(self.outBytes());
+        self.base = try self.rig.scratch(self.outBytes());
+        return true;
+    }
+
+    pub fn run(self: *Pools) !void {
+        const x = self.shape;
+        const t = glm.index.tile(x.n);
+        var a: cuda.Args = .{};
+        inline for (.{ self.in, self.out, @as(i64, @intCast(x.n)), self.pos, @as(c_int, @intCast(x.k)), @as(c_int, @intCast(t.cap)), @as(c_int, @intCast(t.target)), @as(u64, 0) }) |v| a.add(v);
+        try cuda.launch.launch(self.fast, .{ .grid = .{ .x = @intCast(x.rows) }, .block = .{ .x = glm.index.pool_threads } }, self.rig.s, &a);
+        try self.rig.s.synchronize();
+    }
+
+    /// tf_glm_select_rows on the same scores, a 128-thread block a row as cuda_sparse launches it.
+    pub fn replaced(self: *Pools) ![]const u8 {
+        const x = self.shape;
+        var a: cuda.Args = .{};
+        inline for (.{ self.in, self.base, @as(i64, @intCast(x.n)), self.pos, @as(c_int, @intCast(x.k)) }) |v| a.add(v);
+        try cuda.launch.launch(self.narrow, .{ .grid = .{ .x = @intCast(x.rows) }, .block = .{ .x = 128 } }, self.rig.s, &a);
+        return self.rig.bytes(self.base, self.outBytes());
+    }
+
+    pub fn count(self: *Pools) usize {
+        return self.outBytes();
+    }
+
+    pub fn reference(_: *Pools, want: []f64, mag: []f64) !void {
+        @memset(want, 0);
+        @memset(mag, 0);
+    }
+
+    /// tf_glm_select_rows' bytes, then every row's picks against the host's rule over its first max(visible, K) pools.
+    pub fn check(self: *Pools, _: g.Gate) !?[]const u8 {
+        const a = self.rig.gpu.gpa;
+        const x = self.shape;
+        const got = try a.dupe(u8, try self.rig.outputBytes());
+        defer a.free(got);
+        if (std.mem.indexOfDiff(u8, got, try self.replaced())) |at| return try std.fmt.bufPrint(&self.msg, "byte {d} of {d} differs from tf_glm_select_rows'", .{ at, got.len });
+        const keys = try a.alloc(u32, x.n);
+        defer a.free(keys);
+        const idx = try a.alloc(u32, x.n);
+        defer a.free(idx);
+        for (0..x.rows) |r| {
+            const limit = @min(x.n, @max(self.visible(r), x.k));
+            for (keys[0..limit], self.host[r * x.n ..][0..limit]) |*k, s| k.* = orderKey(s);
+            for (idx[0..limit], 0..) |*v, p| v.* = @intCast(p);
+            std.mem.sort(u32, idx[0..limit], keys[0..limit], before);
+            const picks = idx[0..x.k];
+            std.mem.sort(u32, picks, {}, std.sort.asc(u32));
+            for (picks, 0..) |want, i| {
+                const have = std.mem.readInt(i64, got[(r * x.k + i) * 8 ..][0..8], .little);
+                if (have != want) return try std.fmt.bufPrint(&self.msg, "row {d} pick {d}: pool {d}, the host's rule gives {d}", .{ r, i, have, want });
+            }
+        }
+        return null;
+    }
+};
+
+/// The prompt indexer's two kernels, registered together so the runner keeps one entry for them.
+pub const Fast = struct {
+    mod: cuda.Module,
+    scores: @import("contract_glm_index.zig").Scores,
+    pools: Pools,
+
+    /// On the heap, so the kernels' contexts keep their address; registered into `kernels`.
+    pub fn create(a: std.mem.Allocator, gpu: @import("check.zig").Gpu, s: cuda.Stream, fns: glm.glue.Kernels, kernels: *std.ArrayList(g.Kernel)) !*Fast {
+        const idx = @import("contract_glm_index.zig");
+        const self = try a.create(Fast);
+        errdefer a.destroy(self);
+        var m = try cuda.Module.load(gpu.d, cuda.kernels.glm_index);
+        errdefer m.unload();
+        self.* = .{ .mod = m, .scores = idx.Scores{ .rig = .{ .gpu = gpu, .s = s }, .f = fns.scores, .fast = try m.function("tf_glm_index_scores") }, .pools = .{ .rig = .{ .gpu = gpu, .s = s }, .fast = try m.function("tf_glm_index_pools"), .narrow = fns.select_rows } };
+        try kernels.append(a, self.scores.kernel());
+        try kernels.append(a, self.pools.kernel());
+        return self;
+    }
+
+    pub fn destroy(self: *Fast, a: std.mem.Allocator) void {
+        self.pools.deinit();
+        self.scores.deinit();
+        self.mod.unload();
+        a.destroy(self);
+    }
+};
diff --git a/zig/tests/cuda/contract_limits.zig b/zig/tests/cuda/contract_limits.zig
index d8d01b8d74..a9c9690c78 100644
--- a/zig/tests/cuda/contract_limits.zig
+++ b/zig/tests/cuda/contract_limits.zig
@@ -25,7 +25,7 @@ pub fn image(name: []const u8) error{UnknownContractImage}![]const u8 {
     if (std.mem.eql(u8, name, "h3")) return &h3_image.bytes;
     inline for (.{ "qwen4_projection", "qwen4_hc", "qwen4_gdn", "qwen4_ple" }) |n|
         if (std.mem.eql(u8, name, n)) return @field(kk, n);
-    inline for (.{ "tree_round", "act_rows", "gdn_prep", "qk_rows", "norms", "moe", "experts", "qwen36_mtp", "nemotron_attention", "nemotron_route", "nemotron_latent", "qmm_group", "lane_gemv", "prompt16", "qmmf", "ck_experts", "nvfp4_experts", "nvfp4_experts_tile", "rank_sum", "prefill_attention", "vit", "embed_rows", "route_rows", "flashnext_gdn", "flashnext_batch", "flashnext_hc", "flashnext_attn", "flashnext_mtp", "qwen_vision", "exl3_linear", "exl3_experts", "exl3_candidate", "experts3", "experts3_prefill", "drama3", "drama3_fp8", "drama3_codec", "glm_mm", "glm_dflash", "glm_glue", "glm_ops", "nemotron_mamba", "sample", "qwen35", "qwen35_bf16_rows", "qwen35_keep", "qwen35_kv8", "w4a8_prefill", "qwen35_draft", "ck_lane", "gdn_train", "train_q4", "train_attn", "flashnext_train", "train", "attn_train", "qwen27_train", "tree_attention", "block_attention", "window_attention", "window_pages", "moe_glue", "dense_rows", "qmm", "flashnext_fact", "flashnext_hcf", "flashnext_fexp", "experts_rows" }) |n| if (std.mem.eql(u8, name, n)) return @field(kk, n);
+    inline for (.{ "tree_round", "act_rows", "gdn_prep", "qk_rows", "norms", "moe", "experts", "qwen36_mtp", "nemotron_attention", "nemotron_route", "nemotron_latent", "qmm_group", "lane_gemv", "prompt16", "qmmf", "ck_experts", "nvfp4_experts", "nvfp4_experts_tile", "rank_sum", "prefill_attention", "vit", "embed_rows", "route_rows", "flashnext_gdn", "flashnext_batch", "flashnext_hc", "flashnext_attn", "flashnext_mtp", "qwen_vision", "exl3_linear", "exl3_experts", "exl3_candidate", "experts3", "experts3_prefill", "drama3", "drama3_fp8", "drama3_codec", "glm_mm", "glm_dflash", "glm_glue", "glm_ops", "glm_index", "nemotron_mamba", "sample", "qwen35", "qwen35_bf16_rows", "qwen35_keep", "qwen35_kv8", "w4a8_prefill", "qwen35_draft", "ck_lane", "gdn_train", "train_q4", "train_attn", "flashnext_train", "train", "attn_train", "qwen27_train", "tree_attention", "block_attention", "window_attention", "window_pages", "moe_glue", "dense_rows", "qmm", "flashnext_fact", "flashnext_hcf", "flashnext_fexp", "experts_rows" }) |n| if (std.mem.eql(u8, name, n)) return @field(kk, n);
     if (std.mem.eql(u8, name, "contract_controls")) return if (options.controls) &controls.bytes else &.{};
     std.debug.print("FAIL kernel-contract: unknown image name: {s}\n", .{name});
     return error.UnknownContractImage;
diff --git a/zig/tests/cuda/contract_tests.zig b/zig/tests/cuda/contract_tests.zig
index 25e4121ac0..0f6681b598 100644
--- a/zig/tests/cuda/contract_tests.zig
+++ b/zig/tests/cuda/contract_tests.zig
@@ -33,7 +33,6 @@ const latent_k = @import("contract_latent.zig");
 const vit_k = @import("contract_vit_glue.zig");
 const act_round_k = @import("contract_act_round.zig");
 const f32w_k = @import("contract_glm_f32w.zig");
-
 const w4a8_k = @import("contract_w4a8.zig");
 
 const round_k = @import("contract_tree_round.zig");
@@ -538,6 +537,8 @@ pub fn run(gpu: check.Gpu) !void {
     var glm_scores = glm_index_k.Scores{ .rig = .{ .gpu = gpu, .s = s }, .f = glm_fns.scores };
     defer glm_scores.deinit();
     try kernels.append(a, glm_scores.kernel());
+    const glm_fast = try @import("contract_glm_pools.zig").Fast.create(a, gpu, s, glm_fns, &kernels);
+    defer glm_fast.destroy(a);
     var glm_ops = try cuda.Module.load(gpu.d, kk.glm_ops);
     defer glm_ops.unload();
     var glm_tokens = glm_index_k.Tokens{ .rig = .{ .gpu = gpu, .s = s }, .f = try glm_ops.function("tf_glm_dsa_tokens") };
diff --git a/zig/tests/cuda/glm_index_check.zig b/zig/tests/cuda/glm_index_check.zig
new file mode 100644
index 0000000000..48dd2d0ad2
--- /dev/null
+++ b/zig/tests/cuda/glm_index_check.zig
@@ -0,0 +1,269 @@
+//! glm-index-check (GPU): selectTokens with glm_index against tf_glm_scores and tf_glm_select_rows, byte for byte.
+const std = @import("std");
+const cuda = @import("cuda");
+const check = @import("check.zig");
+const glm = @import("glm_cuda");
+const sp = glm.sparse;
+
+const d = 128;
+const page_bytes = 65536; // a page's stride: 128 pools of keys in its first half, as a page that also holds other rows
+
+const Case = struct { cap: usize, pos: usize, rows: usize, heads: usize = 32, paged: bool = false, bits: bool = false };
+const cases = [_]Case{
+    .{ .cap = 4096, .pos = 300, .rows = 512 }, // few visible pools: the rows select over their first 512
+    .{ .cap = 1024, .pos = 4 * 1024 - 520, .rows = 512 },
+    .{ .cap = 4096, .pos = 4 * 4096 - 520, .rows = 512, .bits = true }, // NaN and signed zeros in tied scores
+    .{ .cap = 16384, .pos = 4 * 16384 - 520, .rows = 512 },
+    .{ .cap = 16384, .pos = 4 * 16384 - 1130, .rows = 1100 }, // three score blocks, the later two from s.at
+    .{ .cap = 16384, .pos = 4 * 16384 - 600, .rows = 512, .paged = true },
+    .{ .cap = 4096, .pos = 4 * 4096 - 520, .rows = 512, .heads = 20 }, // fewer heads than a warp's 32 rows of q
+    .{ .cap = 65536, .pos = 4 * 65536 - 520, .rows = 512 },
+    .{ .cap = 147456, .pos = 4 * 147456 - 520, .rows = 512 }, // a cap off the power of two
+    .{ .cap = 262144, .pos = 4 * 262144 - 520, .rows = 512 },
+};
+
+fn bf16(x: f32) u16 {
+    const b: u32 = @bitCast(x);
+    return @intCast((b + 0x7FFF + ((b >> 16) & 1)) >> 16);
+}
+
+const Bufs = std.ArrayList(cuda.DeviceBuffer);
+
+/// n zeroed device bytes, freed with the case.
+fn zeros(gpu: check.Gpu, bufs: *Bufs, n: usize) !u64 {
+    var b = try cuda.DeviceBuffer.alloc(gpu.d, n);
+    errdefer b.free();
+    try b.fill8(0, null);
+    try gpu.d.synchronize(); // before any launch on the non-blocking stream
+    try bufs.append(gpu.gpa, b);
+    return b.ptr;
+}
+
+fn fromHost(gpu: check.Gpu, bufs: *Bufs, bytes: []const u8) !u64 {
+    var b = try cuda.DeviceBuffer.fromHost(gpu.d, bytes);
+    errdefer b.free();
+    try bufs.append(gpu.gpa, b);
+    return b.ptr;
+}
+
+fn scratch(gpu: check.Gpu, bufs: *Bufs, rows: usize, np: usize) !sp.Scratch {
+    const blk: usize = @min(rows, sp.select_rows);
+    return .{ .scores = try zeros(gpu, bufs, blk * np * 4), .ls = try zeros(gpu, bufs, rows * sp.select_list * 4), .li = try zeros(gpu, bufs, rows * sp.select_list * 4), .pools = try zeros(gpu, bufs, rows * sp.topk_pools * 8), .at = try zeros(gpu, bufs, 4) };
+}
+
+fn ms(a: std.Io.Timestamp, b: std.Io.Timestamp) f64 {
+    return @as(f64, @floatFromInt(a.durationTo(b).nanoseconds)) / 1e6;
+}
+
+/// One line before each step, so a failure names the case and step it died in.
+fn step(c: Case, np: usize, what: []const u8) void {
+    std.debug.print("  [np {d} rows {d} heads {d}{s}] {s}\n", .{ np, c.rows, c.heads, if (c.paged) " paged" else "", what });
+}
+
+/// How many rows took each path of glm_index_pools (whole row, candidates, whole row after the candidates missed).
+var paths = [3]usize{ 0, 0, 0 };
+
+pub fn run(gpu: check.Gpu) !void {
+    const gpa = gpu.gpa;
+    var s = try cuda.Stream.init(gpu.d, true);
+    defer s.deinit();
+    var glue_mod = try cuda.Module.load(gpu.d, cuda.kernels.glm_glue);
+    defer glue_mod.unload();
+    var ops = try cuda.Module.load(gpu.d, cuda.kernels.glm_ops);
+    defer ops.unload();
+    const ks = try glm.glue.Kernels.load(glue_mod);
+    var idx = try glm.index.Kernels.load(gpu.d, cuda.kernels.glm_index, .both);
+    defer idx.unload();
+    var idx_scores = idx; // TF_GLM_INDEX_CUDA=scores: the kernels' scores, tf_glm_select_rows' pools
+    idx_scores.select = false;
+    const old: glm.glue.Set = .{ .glue = &ks };
+    const new: glm.glue.Set = .{ .glue = &ks, .index = &idx };
+    const half: glm.glue.Set = .{ .glue = &ks, .index = &idx_scores };
+    var rng = std.Random.DefaultPrng.init(0x105);
+    const rand = rng.random();
+    var fails: usize = 0;
+    for (cases) |c| {
+        const np = sp.poolBucket(c.pos, c.rows, c.cap);
+        const qi = try gpa.alloc(u16, c.rows * c.heads * d);
+        defer gpa.free(qi);
+        for (qi) |*x| x.* = bf16(rand.floatNorm(f32));
+        const wts = try gpa.alloc(u16, c.rows * c.heads);
+        defer gpa.free(wts);
+        for (wts) |*x| x.* = bf16(rand.floatNorm(f32));
+        const pk = try gpa.alloc(u16, c.cap * d);
+        defer gpa.free(pk);
+        for (pk) |*x| x.* = bf16(rand.floatNorm(f32) * 2);
+        var bufs: Bufs = .empty;
+        defer {
+            for (bufs.items) |*b| b.free();
+            bufs.deinit(gpa);
+        }
+        step(c, np, "upload inputs");
+        const dq = try fromHost(gpu, &bufs, std.mem.sliceAsBytes(qi));
+        const dw = try fromHost(gpu, &bufs, std.mem.sliceAsBytes(wts));
+        var table: glm.latent.Table = .{};
+        var dk: u64 = undefined;
+        if (c.paged) {
+            // the keys on shuffled pages, the table long enough for the token lists' latent rows too
+            const page_rows = page_bytes / (512 * 2);
+            const entries = @max((c.cap + 127) / 128, (4 * np + page_rows - 1) / page_rows + 8);
+            const ids = try gpa.alloc(i32, entries);
+            defer gpa.free(ids);
+            for (ids, 0..) |*v, i| v.* = @intCast(i);
+            rand.shuffle(i32, ids);
+            const arena = try gpa.alloc(u8, entries * page_bytes);
+            defer gpa.free(arena);
+            @memset(arena, 0);
+            for (0..c.cap) |p| @memcpy(arena[@as(usize, @intCast(ids[p / 128])) * page_bytes + (p % 128) * d * 2 ..][0 .. d * 2], std.mem.sliceAsBytes(pk[p * d ..][0..d]));
+            dk = try fromHost(gpu, &bufs, arena);
+            table = .{ .pt = try fromHost(gpu, &bufs, std.mem.sliceAsBytes(ids)), .ps = page_bytes, .page_rows = page_rows };
+        } else dk = try fromHost(gpu, &bufs, std.mem.sliceAsBytes(pk));
+        const p0: i32 = @intCast(c.pos);
+        const dpos = try fromHost(gpu, &bufs, std.mem.asBytes(&p0));
+        step(c, np, "allocate scratch");
+        const sa = try scratch(gpu, &bufs, c.rows, np);
+        const sb = try scratch(gpu, &bufs, c.rows, np);
+        const sc = try scratch(gpu, &bufs, c.rows, np);
+        const tok_bytes = c.rows * sp.width * 4;
+        const ta = try zeros(gpu, &bufs, tok_bytes);
+        const tb = try zeros(gpu, &bufs, tok_bytes);
+        const tc = try zeros(gpu, &bufs, tok_bytes);
+        const ca = try zeros(gpu, &bufs, c.rows * 4);
+        const cb = try zeros(gpu, &bufs, c.rows * 4);
+        const cc = try zeros(gpu, &bufs, c.rows * 4);
+        step(c, np, "selectTokens old and new, timed");
+        var ms_old: f64 = 0;
+        var ms_new: f64 = 0;
+        for (0..3) |rep| {
+            const t0 = std.Io.Timestamp.now(gpu.io, .awake);
+            try sp.selectTokens(old, ops, s, dq, dw, c.heads, dk, table, dpos, c.pos, c.rows, c.heads, d, c.cap, sa, ta, ca);
+            try s.synchronize();
+            const t1 = std.Io.Timestamp.now(gpu.io, .awake);
+            try sp.selectTokens(new, ops, s, dq, dw, c.heads, dk, table, dpos, c.pos, c.rows, c.heads, d, c.cap, sb, tb, cb);
+            try s.synchronize();
+            const t2 = std.Io.Timestamp.now(gpu.io, .awake);
+            if (rep > 0) {
+                ms_old += ms(t0, t1) / 2;
+                ms_new += ms(t1, t2) / 2;
+            }
+        }
+        step(c, np, "selectTokens scores-only");
+        try sp.selectTokens(half, ops, s, dq, dw, c.heads, dk, table, dpos, c.pos, c.rows, c.heads, d, c.cap, sc, tc, cc);
+        try s.synchronize();
+        step(c, np, "compare scores, pools, tokens, scores-only");
+        const blk: usize = @min(c.rows, sp.select_rows);
+        fails += try same(gpu, "scores", sa.scores, sb.scores, blk * np * 4, np, c);
+        fails += try same(gpu, "pools", sa.pools, sb.pools, c.rows * sp.topk_pools * 8, np, c);
+        fails += try same(gpu, "tokens", ta, tb, tok_bytes, np, c);
+        fails += try same(gpu, "scores-only scores", sa.scores, sc.scores, blk * np * 4, np, c);
+        fails += try same(gpu, "scores-only pools", sa.pools, sc.pools, c.rows * sp.topk_pools * 8, np, c);
+        fails += try same(gpu, "scores-only tokens", ta, tc, tok_bytes, np, c);
+        std.debug.print("np {d} rows {d} heads {d}{s}: selectTokens {d:.2} ms, with glm_index {d:.2} ms\n", .{ np, c.rows, c.heads, if (c.paged) " paged" else "", ms_old, ms_new });
+        step(c, np, "every score block");
+        fails += try blocks(gpu, &bufs, old, &idx, s, dq, dw, dk, table, c, np);
+        step(c, np, "tied, random and special scores through both selects");
+        fails += try ties(gpu, &bufs, rand, old, &idx, s, sa, sb, dpos, np, c);
+    }
+    std.debug.print("glm_index_pools paths: whole row {d}, candidates {d}, whole row after the candidates missed {d}\n", .{ paths[0], paths[1], paths[2] });
+    try check.expect(paths[0] > 0 and paths[1] > 0 and paths[2] > 0, "glm-index-check: the whole-row, candidate and retry paths must all occur", .{});
+    try check.expect(fails == 0, "glm-index-check: {d} buffers differ", .{fails});
+    check.pass("glm-index-check: scores, pools and token lists byte-equal at {d} cases, tied scores through both selects", .{cases.len});
+}
+
+/// Every 512-row block's scores from tf_glm_scores and from glm_index_scores, as selectTokens launches them.
+fn blocks(gpu: check.Gpu, bufs: *Bufs, set: glm.glue.Set, idx: *const glm.index.Kernels, s: cuda.Stream, dq: u64, dw: u64, dk: u64, t: glm.latent.Table, c: Case, np: usize) !usize {
+    const blk = sp.select_rows;
+    const scale: f32 = @floatCast(std.math.pow(f64, d, -0.5));
+    const loop: usize = if (c.rows >= sp.score_loop_from) sp.score_loop else 1;
+    const first: usize = @min(c.rows, blk);
+    const out_a = try zeros(gpu, bufs, first * np * 4);
+    const out_b = try zeros(gpu, bufs, first * np * 4);
+    var fails: usize = 0;
+    var a: usize = 0;
+    while (a < c.rows) : (a += blk) {
+        const n: usize = @min(blk, c.rows - a);
+        step(c, np, "block: upload its position, clear outputs");
+        const at_v: i32 = @intCast(c.pos + a);
+        const at = try fromHost(gpu, bufs, std.mem.asBytes(&at_v));
+        const nb = (sp.scoreColumns(c.pos + a, n, np) + 63) / 64;
+        try (cuda.DeviceBuffer{ .d = gpu.d, .ptr = out_a, .len = n * np * 4, .borrowed = true }).fill8(0, null);
+        try (cuda.DeviceBuffer{ .d = gpu.d, .ptr = out_b, .len = n * np * 4, .borrowed = true }).fill8(0, null);
+        try gpu.d.synchronize(); // the memsets run on the default stream, not on s
+        step(c, np, "block: tf_glm_scores");
+        var sa: cuda.Args = .{};
+        inline for (.{ dq + a * c.heads * d * 2, dw + a * c.heads * 2, @as(i64, @intCast(c.heads)), dk, out_a, at, @as(i64, @intCast(np)), @as(i64, @intCast(nb)), scale, sp.headScale(c.heads), @as(i64, @intCast(c.heads)), @as(c_int, @intCast(loop)), t.pt, @as(i64, @intCast(t.ps)) }) |v| sa.add(v);
+        try cuda.launch.launch(set.glue.scores, .{ .grid = .{ .x = @intCast(n), .y = @intCast((nb + loop - 1) / loop) }, .block = .{ .x = 128 } }, s, &sa);
+        step(c, np, "block: tf_glm_index_scores");
+        try idx.scores(s, dq + a * c.heads * d * 2, dw + a * c.heads * 2, c.heads, dk, out_b, at, n, np, nb, scale, sp.headScale(c.heads), c.heads, t.pt, t.ps);
+        try s.synchronize();
+        step(c, np, "block: compare");
+        fails += try same(gpu, "block scores", out_a, out_b, n * np * 4, np, c);
+    }
+    return fails;
+}
+
+/// Two device buffers of n bytes compared; 1 if any 4-byte word differs.
+fn same(gpu: check.Gpu, what: []const u8, a: u64, b: u64, n: usize, np: usize, c: Case) !usize {
+    const ha = try gpu.gpa.alloc(u8, n);
+    defer gpu.gpa.free(ha);
+    const hb = try gpu.gpa.alloc(u8, n);
+    defer gpu.gpa.free(hb);
+    std.debug.print("    {s}: read {d} bytes twice\n", .{ what, n });
+    try (cuda.DeviceBuffer{ .d = gpu.d, .ptr = a, .len = n, .borrowed = true }).download(0, ha);
+    try (cuda.DeviceBuffer{ .d = gpu.d, .ptr = b, .len = n, .borrowed = true }).download(0, hb);
+    var diff: usize = 0;
+    for (std.mem.bytesAsSlice(u32, ha), std.mem.bytesAsSlice(u32, hb)) |x, y| diff += @intFromBool(x != y);
+    if (diff != 0) {
+        std.debug.print("FAIL np {d} rows {d} {s}: {d} words differ\n", .{ np, c.rows, what, diff });
+        var shown: usize = 0;
+        for (std.mem.bytesAsSlice(u32, ha), std.mem.bytesAsSlice(u32, hb), 0..) |x, y, i| {
+            if (x == y) continue;
+            std.debug.print("    word {d} (row {d}, column {d} if np wide): {x:0>8} vs {x:0>8}\n", .{ i, i / np, i % np, x, y });
+            shown += 1;
+            if (shown == 6) break;
+        }
+    }
+    return @intFromBool(diff != 0);
+}
+
+/// Tied and random scores (c.bits: with NaN, +-0, +-inf) through both selects, counting glm_index_pools' paths.
+fn ties(gpu: check.Gpu, bufs: *Bufs, rand: std.Random, set: glm.glue.Set, idx: *const glm.index.Kernels, s: cuda.Stream, sa: sp.Scratch, sb: sp.Scratch, at: u64, np: usize, c: Case) !usize {
+    const rows: usize = @min(c.rows, sp.select_rows);
+    const host = try gpu.gpa.alloc(f32, rows * np);
+    defer gpu.gpa.free(host);
+    const path = try zeros(gpu, bufs, rows * 4);
+    const special = [_]f32{ std.math.nan(f32), -std.math.nan(f32), 0.0, -0.0, std.math.inf(f32), -std.math.inf(f32) };
+    var fails: usize = 0;
+    for (0..2) |round| {
+        for (0..rows) |r| {
+            const vis = (c.pos + r + 1) / 4;
+            for (host[r * np ..][0..np], 0..) |*x, i| {
+                x.* = if (i >= vis) -std.math.inf(f32) else if (round == 0) @as(f32, @floatFromInt(rand.uintLessThan(u32, 16))) * 0.25 else rand.floatNorm(f32);
+                if (c.bits and i < vis and rand.uintLessThan(u32, 8) == 0) x.* = special[rand.uintLessThan(usize, special.len)];
+            }
+        }
+        step(c, np, "ties: upload scores");
+        try (cuda.DeviceBuffer{ .d = gpu.d, .ptr = sa.scores, .len = rows * np * 4, .borrowed = true }).upload(0, std.mem.sliceAsBytes(host));
+        step(c, np, "ties: tf_glm_select_rows");
+        var sel: cuda.Args = .{};
+        inline for (.{ sa.scores, sa.pools, @as(i64, @intCast(np)), at, @as(c_int, sp.topk_pools) }) |v| sel.add(v);
+        try cuda.launch.launch(set.glue.select_rows, .{ .grid = .{ .x = @intCast(rows) }, .block = .{ .x = 128 } }, s, &sel);
+        step(c, np, "ties: tf_glm_index_pools");
+        try idx.pools(s, sa.scores, sb.pools, at, rows, np, path);
+        try s.synchronize();
+        step(c, np, "ties: read the paths");
+        var taken: [512]i32 = undefined;
+        try (cuda.DeviceBuffer{ .d = gpu.d, .ptr = path, .len = rows * 4, .borrowed = true }).download(0, std.mem.sliceAsBytes(taken[0..rows]));
+        step(c, np, "ties: count the paths");
+        for (taken[0..rows], 0..) |p, r| {
+            if (p < 0 or p > 2) {
+                std.debug.print("FAIL np {d} row {d}: glm_index_pools reported path {d}\n", .{ np, r, p });
+                return error.TestFailed;
+            }
+            paths[@intCast(p)] += 1;
+        }
+        step(c, np, "ties: compare pools");
+        fails += try same(gpu, if (round == 0) "tied pools" else "random pools", sa.pools, sb.pools, rows * sp.topk_pools * 8, np, c);
+    }
+    return fails;
+}
diff --git a/zig/tests/cuda/main.zig b/zig/tests/cuda/main.zig
index 52c547d425..147128a05a 100644
--- a/zig/tests/cuda/main.zig
+++ b/zig/tests/cuda/main.zig
@@ -35,6 +35,7 @@ const glm_exl3_tests = @import("glm_exl3_tests.zig");
 const glm_vision_tests = @import("glm_vision_tests.zig");
 const glm_mm_fixed = @import("glm_mm_fixed_tests.zig");
 const glm_mm_bench = @import("glm_mm_bench.zig");
+const glm_index_check = @import("glm_index_check.zig");
 const carveout_tests = @import("carveout_tests.zig");
 const nvfp4_tests = @import("nvfp4_tests.zig");
 const modelopt_tests = @import("modelopt_tests.zig");
@@ -228,6 +229,7 @@ const usage =
     \\  glm-vision MODEL F32 BF16 IMG  GLM's image tower and preparation against vision_ref.py's torch runs
     \\  glm-mm-fixed              Prompt K policies: FP64, every row alone/batched, offsets and guarded tails
     \\  glm-mm-bench              GLM's bf16 linears: b16.matmul against glm_mm, by rows (one GPU)
+    \\  glm-index-check           GLM's indexer with glm_index against tf_glm_scores and tf_glm_select_rows, byte for byte (one GPU)
     \\  glm-ref                   GLM projections (matmuls, router, attention, indexer, hc) against float64 (synthetic)
     \\  glm-lanes-cover <ckpt> <dc.json> <rank> <world> <idfile>  copy drafts, cancels, discards, eos and max_new on one lanes.Engine across requests
     \\  glm-lanes-concurrent <ckpt> <dc.json> <rank> <world> <idfile>  concurrent streams on one lanes.Engine, each against its solo run
@@ -364,6 +366,7 @@ fn run(gpu: check.Gpu, cmd: []const u8, rest: []const [:0]const u8) !void {
     if (std.mem.eql(u8, cmd, "glm-exl3")) return glm_exl3_tests.run_(gpu);
     if (std.mem.eql(u8, cmd, "glm-mm-fixed")) return glm_mm_fixed.run(gpu);
     if (std.mem.eql(u8, cmd, "glm-mm-bench")) return glm_mm_bench.run(gpu);
+    if (std.mem.eql(u8, cmd, "glm-index-check")) return glm_index_check.run(gpu);
     if (std.mem.eql(u8, cmd, "glm-vision")) return glm_vision_tests.run(gpu, try arg(rest, 0), try arg(rest, 1), try arg(rest, 2), try arg(rest, 3));
     if (std.mem.eql(u8, cmd, "modelopt")) return modelopt_tests.run(gpu);
     if (std.mem.eql(u8, cmd, "moe")) return moe_tests.run(gpu);