diff --git a/CMakeLists.txt b/CMakeLists.txt
index e63a444..6e1477d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -166,6 +166,11 @@ add_subdirectory(cuBQL)
 
 if (NOT CUBQL_IS_SUBPROJECT)
   add_subdirectory(samples)
+  option(CUBQL_BUILD_TESTS "Build cuBQL correctness tests" OFF)
+  if (CUBQL_BUILD_TESTS)
+    enable_testing()
+    add_subdirectory(tests)
+  endif()
 endif()
 
 #add_subdirectory(testing)
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
new file mode 100644
index 0000000..729d8c1
--- /dev/null
+++ b/tests/CMakeLists.txt
@@ -0,0 +1,24 @@
+# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+# SPDX-License-Identifier: Apache-2.0
+
+if (CUBQL_HAVE_CUDA)
+  add_executable(cubql_warp_unified warp_unified.cu)
+  target_link_libraries(cubql_warp_unified PRIVATE cuBQL)
+  target_compile_definitions(cubql_warp_unified PRIVATE CUBQL_HAVE_CUDA=1)
+  set_target_properties(cubql_warp_unified PROPERTIES
+    CUDA_STANDARD 17 CUDA_STANDARD_REQUIRED ON)
+  add_test(NAME cubql_warp_unified
+    COMMAND cubql_warp_unified --tiles 8 4)
+  add_test(NAME cubql_warp_unified_tails
+    COMMAND cubql_warp_unified --tiles 9 7 --queries 31 15 --window 128)
+  add_test(NAME cubql_warp_unified_small_windows
+    COMMAND cubql_warp_unified --tiles 8 4 --queries 147 74 --stride 7 --window 7)
+  add_test(NAME cubql_warp_unified_regions
+    COMMAND cubql_warp_unified --tiles 64 64 --queries 7 7 --stride 1024 --window 128)
+  add_test(NAME cubql_warp_unified_single_query
+    COMMAND cubql_warp_unified --tiles 1 1 --queries 1 1 --window 512)
+  add_test(NAME cubql_warp_unified_boundaries
+    COMMAND cubql_warp_unified --tiles 3 2 --queries 7 5 --stride 37 --window 1)
+  add_test(NAME cubql_warp_unified_dense
+    COMMAND cubql_warp_unified --tiles 16 16 --queries 15 13 --stride 127 --window 512)
+endif()
diff --git a/tests/data/eda_layout.h b/tests/data/eda_layout.h
new file mode 100644
index 0000000..5ce42e7
--- /dev/null
+++ b/tests/data/eda_layout.h
@@ -0,0 +1,231 @@
+// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+// SPDX-License-Identifier: Apache-2.0
+
+#pragma once
+
+#include "cuBQL/math/box.h"
+#include <algorithm>
+#include <cstddef>
+#include <cstdint>
+#include <fstream>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace eda_layout {
+
+enum class Pattern { Layout, Motif };
+
+struct Rectangle { int x0, y0, x1, y1; };
+
+// Original integer geometry: fingers, routing rails, vertical routes, contacts,
+// and a concentric pair that deliberately produces equal centroid keys.
+inline constexpr Rectangle motif[] = {
+  { 8, 12, 12, 36}, {20, 12, 24, 36}, {32, 12, 36, 36}, {44, 12, 48, 36},
+  {56, 12, 60, 36}, {68, 12, 72, 36}, {80, 12, 84, 36}, {92, 12, 96, 36},
+  { 8, 72, 12, 96}, {20, 72, 24, 96}, {32, 72, 36, 96}, {44, 72, 48, 96},
+  {56, 72, 60, 96}, {68, 72, 72, 96}, {80, 72, 84, 96}, {92, 72, 96, 96},
+  { 4,  4,100,  8}, { 4, 40,100, 44}, { 4, 64,100, 68}, { 4,100,100,104},
+  { 2,  4,  4,104}, {104,  4,106,104}, {16, 44, 18, 64}, {88, 44, 90, 64},
+  { 8, 46, 14, 52}, {26, 46, 32, 52}, {44, 46, 50, 52}, {62, 46, 68, 52},
+  {80, 46, 86, 52}, {98, 46,104, 52}, {48, 54, 56, 62}, {50, 56, 54, 60}
+};
+inline constexpr uint32_t rectanglesPerTile = sizeof(motif)/sizeof(motif[0]);
+inline constexpr int tilePitch = 128;
+
+inline Rectangle layoutRectangle(uint32_t x, uint32_t y, uint32_t r)
+{
+  const uint32_t groupX = x/3;
+  const uint32_t groupY = y/2;
+  const uint32_t selection = (5*groupX+7*groupY+3*(groupX^groupY))%19;
+  uint32_t kind = selection < 8 ? 0 : selection < 13 ? 1 : selection < 18 ? 2 : 3;
+  if ((x%11 == 5 && y%5 != 1) || (y%9 == 4 && x%7 != 2)) kind = 3;
+  Rectangle rect;
+  if (kind == 0) {
+    if (r < 16) {
+      const int x0 = 8+14*int(r%8);
+      const int y0 = 8+72*int(r/8);
+      rect = {x0,y0,x0+10,y0+34};
+    } else if (r < 20) {
+      const int i = int(r-16);
+      const int y0 = 4+72*(i/2)+38*(i%2);
+      rect = {2,y0,122,y0+4};
+    } else if (r < 28) {
+      const int x0 = 10+14*int(r-20);
+      rect = {x0,54,x0+8,62};
+    } else if (r < 30) {
+      const int x0 = 2+118*int(r-28);
+      rect = {x0,4,x0+4,118};
+    } else {
+      rect = r == 30 ? Rectangle{56,64,72,74} : Rectangle{58,66,70,72};
+    }
+  } else if (kind == 1) {
+    if (r < 16) {
+      const int group = int(r/4);
+      const int x0 = 8+60*(group%2)+12*int(r%4);
+      const int y0 = 8+60*(group/2);
+      rect = {x0,y0,x0+8,y0+28};
+    } else if (r < 24) {
+      const int i = int(r-16);
+      const int x0 = 4+60*((i/2)%2);
+      const int y0 = 4+60*(i/4)+32*(i%2);
+      rect = {x0,y0,x0+52,y0+4};
+    } else if (r < 28) {
+      const int i = int(r-24);
+      const int x0 = 4+60*(i%2);
+      const int y0 = 4+60*(i/2);
+      rect = {x0,y0,x0+4,y0+36};
+    } else if (r < 30) {
+      const int y0 = 48+72*int(r-28);
+      rect = {2,y0,124,y0+4};
+    } else {
+      rect = r == 30 ? Rectangle{56,48,72,60} : Rectangle{58,50,70,58};
+    }
+  } else if (kind == 2) {
+    if (r < 16) {
+      const int i = int(r/2);
+      const int x0 = 2+30*(i%4);
+      const int y0 = 8+60*(i/4);
+      rect = r%2 == 0 ? Rectangle{x0,y0,x0+24,y0+10}
+                      : Rectangle{x0+16,y0+10,x0+24,y0+38};
+    } else if (r < 24) {
+      const int i = int(r-16);
+      const int x0 = 8+30*(i%4);
+      const int y0 = 28+60*(i/4);
+      rect = {x0,y0,x0+8,y0+8};
+    } else if (r < 28) {
+      const int i = int(r-24);
+      const int y0 = 4+56*(i/2)+4*(i%2);
+      rect = {4,y0,124,y0+3};
+    } else if (r < 30) {
+      const int x0 = 4+116*int(r-28);
+      rect = {x0,4,x0+4,120};
+    } else {
+      rect = r == 30 ? Rectangle{56,50,72,58} : Rectangle{58,52,70,56};
+    }
+  } else {
+    if (r < 16) {
+      const int i = int(r/2);
+      const int x0 = 8+72*(i%2);
+      const int y0 = 8+28*(i/2);
+      rect = r%2 == 0 ? Rectangle{x0,y0,x0+36,y0+6}
+                      : Rectangle{x0+28,y0+6,x0+36,y0+20};
+    } else if (r < 28) {
+      const int i = int(r-16);
+      const int x0 = 12+12*(i%3)+72*(i/6);
+      const int y0 = 30+32*((i%6)/3);
+      rect = {x0,y0,x0+6,y0+6};
+    } else if (r < 30) {
+      const int x0 = 2+120*int(r-28);
+      rect = {x0,4,x0+4,124};
+    } else {
+      rect = r == 30 ? Rectangle{58,58,70,70} : Rectangle{60,60,68,68};
+    }
+  }
+  if ((groupX+groupY)%2) {
+    const int x0 = tilePitch-rect.x1;
+    rect.x1 = tilePitch-rect.x0;
+    rect.x0 = x0;
+  }
+  if ((x/5+y/3)%2) {
+    const int y0 = tilePitch-rect.y1;
+    rect.y1 = tilePitch-rect.y0;
+    rect.y0 = y0;
+  }
+  return rect;
+}
+
+template<typename T, int D>
+std::vector<cuBQL::box_t<T,D>> makeBoxes(uint32_t tilesX, uint32_t tilesY,
+                                      Pattern pattern=Pattern::Layout)
+{
+  static_assert(D >= 2 && D <= 4);
+  if (tilesX == 0 || tilesY == 0 || tilesX > 131072 || tilesY > 131072)
+    throw std::invalid_argument("Synthetic layout requires 1 to 131072 tiles per axis");
+  if (uint64_t(tilesX)*tilesY*rectanglesPerTile > uint64_t(std::numeric_limits<int>::max())/2)
+    throw std::invalid_argument("Synthetic layout exceeds the radix builder node limit");
+  std::vector<cuBQL::box_t<T,D>> boxes(std::size_t(tilesX)*tilesY*rectanglesPerTile);
+  std::size_t index = 0;
+  for (uint32_t y=0;y<tilesY;++y)
+    for (uint32_t x=0;x<tilesX;++x)
+      for (uint32_t r=0;r<rectanglesPerTile;++r) {
+        const Rectangle rect = pattern == Pattern::Motif ? motif[r] : layoutRectangle(x,y,r);
+        const int dx = int(x)*tilePitch-tilePitch;
+        const int dy = int(y)*tilePitch-tilePitch;
+        auto &box = boxes[index++];
+        box.lower[0] = T(dx+rect.x0); box.upper[0] = T(dx+rect.x1);
+        box.lower[1] = T(dy+rect.y0); box.upper[1] = T(dy+rect.y1);
+        for (int axis=2;axis<D;++axis) {
+          const int lower = axis == 2 ? 4*int(r%4)-8 : 4*int(r/4)-16;
+          box.lower[axis] = T(lower);
+          box.upper[axis] = T(lower+2);
+        }
+      }
+  return boxes;
+}
+
+// Canonical little-endian signed integer coordinates, independent of scalar type.
+template<typename T, int D>
+uint64_t fingerprint(const std::vector<cuBQL::box_t<T,D>> &boxes)
+{
+  uint64_t hash = UINT64_C(14695981039346656037);
+  for (const auto &box : boxes)
+    for (int axis=0;axis<D;++axis)
+      for (int side=0;side<2;++side) {
+        uint64_t value = uint64_t(int64_t(side ? box.upper[axis] : box.lower[axis]));
+        for (int byte=0;byte<8;++byte) {
+          hash ^= (value >> (8*byte)) & 255u;
+          hash *= UINT64_C(1099511628211);
+        }
+      }
+  return hash;
+}
+
+inline void writeSvg(const char *path, uint32_t tilesX=16, uint32_t tilesY=8,
+                     Pattern pattern=Pattern::Layout)
+{
+  if (!path || !*path)
+    throw std::invalid_argument("SVG output path must not be empty");
+  if (tilesX == 0 || tilesY == 0 || tilesX > 32 || tilesY > 32)
+    throw std::invalid_argument("SVG preview requires 1 to 32 tiles per axis");
+  const auto boxes = makeBoxes<float,2>(tilesX,tilesY,pattern);
+  auto bounds = boxes.front();
+  for (const auto &box : boxes)
+    bounds.grow(box);
+  const int left = int(bounds.lower[0])-18;
+  const int top = int(bounds.lower[1])-46;
+  const int width = std::max(340,int(bounds.upper[0]-bounds.lower[0])+36);
+  const int height = int(bounds.upper[1]-bounds.lower[1])+64;
+  std::ofstream out;
+  out.exceptions(std::ios::failbit | std::ios::badbit);
+  try {
+    out.open(path);
+    out << "<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\""
+        << left << ' ' << top << ' ' << width << ' ' << height << "\">\n"
+        << "<title>Original synthetic layout</title>\n"
+        << "<rect x=\"" << left << "\" y=\"" << top
+        << "\" width=\"" << width << "\" height=\"" << height
+        << "\" fill=\"#000000\"/>\n"
+        << "<text x=\"" << int(bounds.lower[0]) << "\" y=\""
+        << int(bounds.lower[1])-28
+        << "\" fill=\"#dbe5eb\" font-family=\"sans-serif\" font-size=\"14\">"
+        << "Original synthetic layout: " << boxes.size() << " rectangles</text>\n";
+    for (std::size_t i=0;i<boxes.size();++i) {
+      const auto &box = boxes[i];
+      const uint32_t r = uint32_t(i%rectanglesPerTile);
+      const char *color = pattern == Pattern::Layout ? "#72b59f"
+                          : (r >= 16 && r < 24) || r == 31 ? "#52c7bd" : "#e7be62";
+      out << "<rect x=\"" << int(box.lower[0]) << "\" y=\"" << int(box.lower[1])
+          << "\" width=\"" << int(box.upper[0]-box.lower[0])
+          << "\" height=\"" << int(box.upper[1]-box.lower[1])
+          << "\" fill=\"" << color << "\"/>\n";
+    }
+    out << "</svg>\n";
+    out.close();
+  } catch (const std::ios_base::failure &) {
+    throw std::runtime_error(std::string("Failed to write SVG preview: ")+path);
+  }
+}
+
+} // namespace eda_layout
diff --git a/tests/data/sample.jpg b/tests/data/sample.jpg
new file mode 100644
index 0000000..16623a5
Binary files /dev/null and b/tests/data/sample.jpg differ
diff --git a/tests/warp_unified.cu b/tests/warp_unified.cu
new file mode 100644
index 0000000..e5507bd
--- /dev/null
+++ b/tests/warp_unified.cu
@@ -0,0 +1,716 @@
+// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+// SPDX-License-Identifier: Apache-2.0
+
+#include <cuda_runtime.h>
+#include "cuBQL/bvh.h"
+#include "cuBQL/builder/cuda.h"
+#include "cuBQL/builder/cuda/radix.h"
+#include "cuBQL/traversal/fixedAnyShapeQuery.h"
+#include "data/eda_layout.h"
+#include "warp_unified.cuh"
+
+#include <algorithm>
+#include <cerrno>
+#include <climits>
+#include <cmath>
+#include <cstdint>
+#include <cstdlib>
+#include <cstring>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace {
+  using warpUnified::Admin;
+  using warpUnified::box2f;
+  using warpUnified::bvh2f;
+  constexpr uint32_t blockSize = 64;
+  constexpr uint32_t referenceBlockSize = 128;
+  constexpr uint32_t referenceStackDepth = 64;
+  constexpr uint32_t exactIntegerLimit = 1u<<24;
+  constexpr int warmupCount = 2;
+  constexpr int roundCount = 12;
+  const char *variantNames[] = {"reference","before","after"};
+
+  void check(cudaError_t error)
+  {
+    if (error != cudaSuccess)
+      throw std::runtime_error(cudaGetErrorString(error));
+  }
+
+  template<typename T>
+  struct DeviceBuffer {
+    T *data = nullptr;
+
+    explicit DeviceBuffer(size_t count)
+    {
+      if (count > std::numeric_limits<size_t>::max()/sizeof(T))
+        throw std::runtime_error("device allocation size overflow");
+      if (count)
+        check(cudaMalloc(reinterpret_cast<void **>(&data),count*sizeof(T)));
+    }
+    ~DeviceBuffer() { cudaFree(data); }
+    DeviceBuffer(const DeviceBuffer &) = delete;
+    DeviceBuffer &operator=(const DeviceBuffer &) = delete;
+  };
+
+  struct Event {
+    cudaEvent_t handle = nullptr;
+    Event() { check(cudaEventCreate(&handle)); }
+    ~Event() { cudaEventDestroy(handle); }
+    Event(const Event &) = delete;
+    Event &operator=(const Event &) = delete;
+  };
+
+  struct Tree {
+    bvh2f bvh;
+    cuBQL::DeviceMemoryResource resource;
+    ~Tree()
+    {
+      cudaFree(bvh.nodes);
+      cudaFree(bvh.primIDs);
+    }
+    Tree() = default;
+    Tree(const Tree &) = delete;
+    Tree &operator=(const Tree &) = delete;
+  };
+
+  struct Options {
+    uint32_t tilesX = 8, tilesY = 4;
+    uint32_t queriesX = 0, queriesY = 0;
+    uint32_t stride = 32, window = 100;
+    bool benchmark = false, help = false;
+  };
+
+  struct Grid {
+    uint32_t x, y, stride, window;
+    int originX = -128, originY = -128;
+  };
+
+  struct QueryResult {
+    uint64_t area = 0;
+    uint32_t hits = 0;
+  };
+
+  uint32_t positiveInteger(const char *text)
+  {
+    errno = 0;
+    char *end = nullptr;
+    const auto value = std::strtoull(text,&end,10);
+    if (!*text || *text == '-' || errno || *end || !value || value > INT_MAX)
+      throw std::runtime_error("expected a positive integer <= INT_MAX");
+    return uint32_t(value);
+  }
+
+  Options parseOptions(int argc, char **argv)
+  {
+    Options options;
+    bool explicitTiles = false;
+    bool explicitStride = false, explicitWindow = false;
+    for (int i=1;i<argc;++i) {
+      const std::string arg = argv[i];
+      if (arg == "--help")
+        options.help = true;
+      else if (arg == "--benchmark")
+        options.benchmark = true;
+      else if ((arg == "--tiles" || arg == "--queries") && i+2 < argc) {
+        const uint32_t x = positiveInteger(argv[++i]);
+        const uint32_t y = positiveInteger(argv[++i]);
+        if (arg == "--tiles") {
+          options.tilesX = x;
+          options.tilesY = y;
+          explicitTiles = true;
+        } else {
+          options.queriesX = x;
+          options.queriesY = y;
+        }
+      } else if (arg == "--stride" && i+1 < argc) {
+        options.stride = positiveInteger(argv[++i]);
+        explicitStride = true;
+      } else if (arg == "--window" && i+1 < argc) {
+        options.window = positiveInteger(argv[++i]);
+        explicitWindow = true;
+      } else
+        throw std::runtime_error("unknown option or missing argument: "+arg);
+    }
+    if (options.help)
+      return options;
+    if (options.benchmark) {
+      if (!explicitTiles) {
+        options.tilesX = 256;
+        options.tilesY = 256;
+      }
+      if (!explicitStride)
+        options.stride = 216;
+      if (!explicitWindow)
+        options.window = 1792;
+    }
+    if (options.stride > exactIntegerLimit || options.window > exactIntegerLimit)
+      throw std::runtime_error("query parameters exceed exact float-coordinate limits");
+    if (options.tilesX > 131072 || options.tilesY > 131072)
+      throw std::runtime_error("layout requires at most 131072 tiles per axis");
+    const uint64_t boxCount = uint64_t(options.tilesX)*options.tilesY*eda_layout::rectanglesPerTile;
+    if (boxCount > uint64_t(INT_MAX)/2
+        || boxCount > std::numeric_limits<size_t>::max()/sizeof(box2f))
+      throw std::runtime_error("layout exceeds builder node or allocation limits");
+    if (!options.queriesX) {
+      const uint64_t extentX = uint64_t(options.tilesX)*128;
+      const uint64_t extentY = uint64_t(options.tilesY)*128;
+      if (options.benchmark) {
+        options.queriesX = uint32_t(1+(extentX > options.window ? extentX-options.window : 0)/options.stride);
+        options.queriesY = uint32_t(1+(extentY > options.window ? extentY-options.window : 0)/options.stride);
+      } else {
+        options.queriesX = uint32_t((extentX+options.stride-1)/options.stride);
+        options.queriesY = uint32_t((extentY+options.stride-1)/options.stride);
+      }
+    }
+    const uint64_t count = uint64_t(options.queriesX)*options.queriesY;
+    if (!count || count > INT_MAX
+        || count > std::numeric_limits<size_t>::max()/sizeof(uint32_t))
+      throw std::runtime_error("query grid exceeds indexing or allocation limits");
+    const int64_t endX = -128+int64_t(options.queriesX-1)*options.stride+options.window;
+    const int64_t endY = -128+int64_t(options.queriesY-1)*options.stride+options.window;
+    if (endX > exactIntegerLimit || endY > exactIntegerLimit)
+      throw std::runtime_error("query endpoints exceed exact float integer range");
+    return options;
+  }
+
+  void printHelp()
+  {
+    std::cout << "Usage: cubql_warp_unified [--benchmark] [--tiles X Y] [--queries X Y]\n"
+                 "                          [--stride N] [--window N]\n"
+                 "Defaults: 8x4 tiles, stride 32, window 100, origin -128,-128.\n"
+                 "--benchmark uses the Morton experiment's 256x256 layout (2097152 rectangles).\n"
+                 "Benchmark defaults: 144x144 overlapping queries, stride 216, window 1792.\n"
+                 "Compares scalar radix leaf1, Andrea WarpUnified, and the lane-0 fix.\n"
+                 "The two WarpUnified variants share a radix leaf16 BVH and Morton query order.\n";
+  }
+
+  __host__ __device__ box2f queryBox(uint32_t id, Grid grid)
+  {
+    const int64_t x = int64_t(grid.originX)+int64_t(id%grid.x)*grid.stride;
+    const int64_t y = int64_t(grid.originY)+int64_t(id/grid.x)*grid.stride;
+    return box2f{cuBQL::vec2f{float(x),float(y)},
+                 cuBQL::vec2f{float(x+grid.window),float(y+grid.window)}};
+  }
+
+  // Scalar baseline from window_test.cu at 396074e.
+  template<typename Lambda>
+  inline __cubql_both void forEachLeaf2D(const Lambda &lambda_to_call_on_each_leaf,
+                                        const bvh2f bvh, const box2f query_box)
+  {
+    Admin traversal_stack[referenceStackDepth], *stack_ptr = traversal_stack;
+    Admin node = bvh.nodes[0].admin;
+    while (true) {
+      while (true) {
+        if (node.count != 0) {
+          break;
+        }
+        const uint32_t n0_idx = static_cast<uint32_t>(node.offset)+0;
+        const uint32_t n1_idx = static_cast<uint32_t>(node.offset)+1;
+        const bvh2f::node_t n0 = bvh.nodes[n0_idx];
+        const bvh2f::node_t n1 = bvh.nodes[n1_idx];
+        const bool o0 = query_box.overlaps(n0.bounds);
+        const bool o1 = query_box.overlaps(n1.bounds);
+        if (o0) {
+          if (o1) {
+            *stack_ptr++ = n1.admin;
+          }
+          node = n0.admin;
+        } else {
+          if (o1) {
+            node = n1.admin;
+          } else {
+            node.count = 0;
+            break;
+          }
+        }
+      }
+      if (node.count != 0) {
+        const int leaf_result = lambda_to_call_on_each_leaf(
+          bvh.primIDs+node.offset,static_cast<uint32_t>(node.count));
+        if (leaf_result == CUBQL_TERMINATE_TRAVERSAL) {
+          return;
+        }
+      }
+      if (stack_ptr == traversal_stack) {
+        return;
+      }
+      node = *--stack_ptr;
+    }
+  }
+
+  template<typename Lambda>
+  inline __cubql_both void forEachPrim2D(const Lambda &lambda_to_call_on_each_prim,
+                                        const bvh2f bvh, const box2f query_box)
+  {
+    auto leaf_code = [&lambda_to_call_on_each_prim]
+      (const uint32_t *prim_ids, size_t num_prims) -> int {
+        for (int i=0;i<static_cast<int>(num_prims);++i) {
+          if (lambda_to_call_on_each_prim(prim_ids[i]) == CUBQL_TERMINATE_TRAVERSAL) {
+            return CUBQL_TERMINATE_TRAVERSAL;
+          }
+        }
+        return CUBQL_CONTINUE_TRAVERSAL;
+      };
+    forEachLeaf2D(leaf_code,bvh,query_box);
+  }
+
+  struct WindowAreaAccumulator {
+    const box2f *boxes = nullptr;
+    box2f window;
+    mutable float total_area = 0.f;
+
+    __device__ int operator()(uint32_t box_id) const
+    {
+      const box2f box = boxes[box_id];
+      const float x0 = fmaxf(window.lower.x,box.lower.x);
+      const float y0 = fmaxf(window.lower.y,box.lower.y);
+      const float x1 = fminf(window.upper.x,box.upper.x);
+      const float y1 = fminf(window.upper.y,box.upper.y);
+      const float dx = fmaxf(0.f,x1-x0);
+      const float dy = fmaxf(0.f,y1-y0);
+      total_area += dx*dy;
+      return CUBQL_CONTINUE_TRAVERSAL;
+    }
+  };
+
+  __global__ void referenceQueryKernel(bvh2f bvh, const box2f *boxes, Grid grid,
+                                      uint32_t count, float *results)
+  {
+    const uint64_t id = uint64_t(blockIdx.x)*blockDim.x+threadIdx.x;
+    if (id >= count)
+      return;
+    if (!bvh.numPrims) {
+      results[id] = 0.f;
+      return;
+    }
+    WindowAreaAccumulator accumulator{boxes,queryBox(uint32_t(id),grid),0.f};
+    forEachPrim2D(accumulator,bvh,accumulator.window);
+    results[id] = accumulator.total_area;
+  }
+
+  template<bool Lane0>
+  __global__ void queryKernel(bvh2f bvh, const box2f *boxes, Grid grid,
+                             const uint32_t *dispatch, uint32_t count, float *results)
+  {
+    const uint64_t warp = (uint64_t(blockIdx.x)*blockDim.x+threadIdx.x)/warpUnified::warpSize;
+    if (warp >= count)
+      return;
+    const uint32_t id = __ldg(&dispatch[warp]);
+    Admin stack[warpUnified::stackDepth];
+    const float area = warpUnified::queryArea<Lane0>(bvh,queryBox(id,grid),boxes,stack);
+    if ((threadIdx.x&(warpUnified::warpSize-1)) == 0)
+      results[id] = area;
+  }
+
+  void appendMorton(std::vector<uint32_t> &ids, uint32_t x, uint32_t y,
+                    uint32_t side, Grid grid)
+  {
+    if (x >= grid.x || y >= grid.y)
+      return;
+    if (side == 1) {
+      ids.push_back(y*grid.x+x);
+      return;
+    }
+    const uint32_t half = side/2;
+    appendMorton(ids,x,y,half,grid);
+    appendMorton(ids,x+half,y,half,grid);
+    appendMorton(ids,x,y+half,half,grid);
+    appendMorton(ids,x+half,y+half,half,grid);
+  }
+
+  std::vector<uint32_t> makeDispatch(Grid grid)
+  {
+    const uint32_t count = grid.x*grid.y;
+    std::vector<uint32_t> ids;
+    ids.reserve(count);
+    uint32_t side = 1;
+    while (side < std::max(grid.x,grid.y))
+      side *= 2;
+    appendMorton(ids,0,0,side,grid);
+    if (ids.size() != count)
+      throw std::runtime_error("Morton dispatch misses queries");
+    std::vector<unsigned char> seen(count,0);
+    for (const uint32_t id : ids)
+      if (id >= count || seen[id]++)
+        throw std::runtime_error("Morton dispatch is not a permutation");
+    return ids;
+  }
+
+  void launch(int variant, bvh2f reference, bvh2f optimized, const box2f *boxes,
+              const box2f *reordered, Grid grid, const uint32_t *dispatch, float *results)
+  {
+    const uint32_t count = grid.x*grid.y;
+    const uint32_t blocks = (count+1)/2;
+    if (variant == 0)
+      referenceQueryKernel<<<uint32_t((uint64_t(count)+referenceBlockSize-1)/referenceBlockSize),referenceBlockSize>>>(
+        reference,boxes,grid,count,results);
+    else if (variant == 1)
+      queryKernel<false><<<blocks,blockSize>>>(optimized,reordered,grid,dispatch,count,results);
+    else
+      queryKernel<true><<<blocks,blockSize>>>(optimized,reordered,grid,dispatch,count,results);
+    check(cudaGetLastError());
+  }
+
+  bool contains(const box2f &outer, const box2f &inner)
+  {
+    return outer.lower.x <= inner.lower.x && outer.lower.y <= inner.lower.y
+      && outer.upper.x >= inner.upper.x && outer.upper.y >= inner.upper.y;
+  }
+
+  bool validBounds(const box2f &box)
+  {
+    return std::isfinite(box.lower.x) && std::isfinite(box.lower.y)
+      && std::isfinite(box.upper.x) && std::isfinite(box.upper.y)
+      && box.lower.x <= box.upper.x && box.lower.y <= box.upper.y;
+  }
+
+  void auditTree(bvh2f bvh, const std::vector<box2f> &boxes, const box2f *reordered,
+                 uint32_t stackCapacity)
+  {
+    if (!bvh.nodes || !bvh.primIDs || !bvh.numNodes || bvh.numPrims != boxes.size())
+      throw std::runtime_error("unexpected BVH counts or storage");
+    std::vector<bvh2f::node_t> nodes(bvh.numNodes);
+    std::vector<uint32_t> ids(bvh.numPrims);
+    std::vector<box2f> gathered(reordered ? boxes.size() : 0);
+    check(cudaMemcpy(nodes.data(),bvh.nodes,nodes.size()*sizeof(nodes[0]),cudaMemcpyDeviceToHost));
+    check(cudaMemcpy(ids.data(),bvh.primIDs,ids.size()*sizeof(ids[0]),cudaMemcpyDeviceToHost));
+    if (reordered)
+      check(cudaMemcpy(gathered.data(),reordered,gathered.size()*sizeof(gathered[0]),cudaMemcpyDeviceToHost));
+    std::vector<unsigned char> primitiveSeen(boxes.size(),0), coverage(boxes.size(),0);
+    for (size_t i=0;i<ids.size();++i) {
+      if (ids[i] >= boxes.size() || primitiveSeen[ids[i]]++)
+        throw std::runtime_error("BVH primitive IDs are not a permutation");
+      if (reordered) {
+        const auto &expected = boxes[ids[i]];
+        const auto &actual = gathered[i];
+        if (actual.lower.x != expected.lower.x || actual.lower.y != expected.lower.y
+            || actual.upper.x != expected.upper.x || actual.upper.y != expected.upper.y)
+          throw std::runtime_error("leaf-order box gather differs from primitive IDs");
+      }
+    }
+    struct Pending { uint32_t index, depth; };
+    std::vector<Pending> pending{{0,0}};
+    std::vector<unsigned char> seen(nodes.size(),0);
+    uint32_t maxDepth = 0;
+    while (!pending.empty()) {
+      const Pending current = pending.back();
+      pending.pop_back();
+      if (current.index >= nodes.size() || seen[current.index]++)
+        throw std::runtime_error("invalid or repeated reachable BVH node");
+      const auto &node = nodes[current.index];
+      if (!validBounds(node.bounds))
+        throw std::runtime_error("invalid BVH node bounds");
+      maxDepth = std::max(maxDepth,current.depth);
+      if (current.depth > stackCapacity)
+        throw std::runtime_error("BVH exceeds traversal stack capacity");
+      const Admin admin = node.admin;
+      if (admin.count) {
+        if (admin.offset+admin.count > bvh.numPrims)
+          throw std::runtime_error("invalid BVH leaf range");
+        for (uint64_t i=admin.offset;i<admin.offset+admin.count;++i)
+          if (coverage[size_t(i)]++
+              || !contains(node.bounds,reordered ? gathered[size_t(i)] : boxes[ids[size_t(i)]]))
+            throw std::runtime_error("overlapping leaf ranges or insufficient leaf bounds");
+      } else {
+        if (admin.offset < 2 || (admin.offset&1) || admin.offset+1 >= nodes.size())
+          throw std::runtime_error("invalid BVH child pair");
+        for (uint32_t child=0;child<2;++child) {
+          const uint32_t index = uint32_t(admin.offset+child);
+          if (!contains(node.bounds,nodes[index].bounds))
+            throw std::runtime_error("BVH parent does not contain child bounds");
+          pending.push_back({index,current.depth+1});
+        }
+      }
+    }
+    if (std::find(coverage.begin(),coverage.end(),0) != coverage.end())
+      throw std::runtime_error("BVH leaves do not cover every primitive");
+    std::cout << "tree_audit=PASS max_depth=" << maxDepth
+              << " stack_capacity=" << stackCapacity << '\n';
+  }
+
+  QueryResult bruteQuery(uint32_t id, const std::vector<box2f> &boxes, Grid grid)
+  {
+    const auto query = queryBox(id,grid);
+    QueryResult result;
+    for (const auto &box : boxes) {
+      const int64_t dx = int64_t(std::min(box.upper.x,query.upper.x))
+        -int64_t(std::max(box.lower.x,query.lower.x));
+      const int64_t dy = int64_t(std::min(box.upper.y,query.upper.y))
+        -int64_t(std::max(box.lower.y,query.lower.y));
+      if (dx > 0 && dy > 0) {
+        result.area += uint64_t(dx*dy);
+        ++result.hits;
+      }
+    }
+    return result;
+  }
+
+  QueryResult queryOracle(int64_t x0, int64_t y0, uint32_t window,
+                          uint32_t tilesX, uint32_t tilesY)
+  {
+    QueryResult result;
+    const int64_t x1 = x0+window, y1 = y0+window;
+    const int64_t lowerX = std::max(x0,int64_t(-128));
+    const int64_t lowerY = std::max(y0,int64_t(-128));
+    const int64_t upperX = std::min(x1,int64_t(tilesX)*eda_layout::tilePitch-128);
+    const int64_t upperY = std::min(y1,int64_t(tilesY)*eda_layout::tilePitch-128);
+    if (upperX <= lowerX || upperY <= lowerY)
+      return result;
+    const uint32_t firstX = uint32_t((lowerX+128)/eda_layout::tilePitch);
+    const uint32_t firstY = uint32_t((lowerY+128)/eda_layout::tilePitch);
+    const uint32_t lastX = uint32_t((upperX+127)/eda_layout::tilePitch);
+    const uint32_t lastY = uint32_t((upperY+127)/eda_layout::tilePitch);
+    for (uint32_t y=firstY;y<=lastY;++y)
+      for (uint32_t x=firstX;x<=lastX;++x)
+        for (uint32_t r=0;r<eda_layout::rectanglesPerTile;++r) {
+          const auto rect = eda_layout::layoutRectangle(x,y,r);
+          const int64_t originX = int64_t(x)*eda_layout::tilePitch-128;
+          const int64_t originY = int64_t(y)*eda_layout::tilePitch-128;
+          const int64_t dx = std::min(originX+rect.x1,x1)-std::max(originX+rect.x0,x0);
+          const int64_t dy = std::min(originY+rect.y1,y1)-std::max(originY+rect.y0,y0);
+          if (dx > 0 && dy > 0) {
+            result.area += uint64_t(dx*dy);
+            ++result.hits;
+          }
+        }
+    return result;
+  }
+
+  std::vector<float> expectedResults(const Options &options,
+                                     const std::vector<box2f> &boxes, Grid grid)
+  {
+    const uint32_t count = grid.x*grid.y;
+    std::vector<float> expected(count);
+    const bool brute = uint64_t(count)*boxes.size() <= UINT64_C(8000000);
+    uint64_t totalHits = 0, emptyQueries = 0;
+    uint32_t maxHits = 0;
+    for (uint32_t id=0;id<count;++id) {
+      const auto result = queryOracle(
+        grid.originX+int64_t(id%grid.x)*grid.stride,
+        grid.originY+int64_t(id/grid.x)*grid.stride,
+        grid.window,options.tilesX,options.tilesY);
+      if (result.area > exactIntegerLimit)
+        throw std::runtime_error("query area exceeds exact float integer range");
+      expected[id] = float(result.area);
+      if (brute) {
+        const auto independent = bruteQuery(id,boxes,grid);
+        if (independent.area != result.area || independent.hits != result.hits)
+          throw std::runtime_error("integer fixture oracle differs from brute force");
+      }
+      totalHits += result.hits;
+      maxHits = std::max(maxHits,result.hits);
+      emptyQueries += result.hits == 0;
+    }
+    std::cout << "cpu_oracle=all_queries brute_force=" << (brute ? "all_queries_PASS" : "small_tests_only")
+              << " empty_queries=" << emptyQueries << " mean_hits=" << double(totalHits)/count
+              << " max_hits=" << maxHits << '\n';
+    return expected;
+  }
+
+  uint64_t outputFingerprint(const std::vector<float> &values)
+  {
+    uint64_t hash = UINT64_C(14695981039346656037);
+    for (const float value : values) {
+      uint32_t bits = 0;
+      std::memcpy(&bits,&value,sizeof(bits));
+      for (int byte=0;byte<4;++byte) {
+        hash ^= (bits>>(8*byte))&255u;
+        hash *= UINT64_C(1099511628211);
+      }
+    }
+    return hash;
+  }
+
+  void validate(const std::vector<float> &expected, const float *device,
+                int variant, const char *phase)
+  {
+    std::vector<float> actual(expected.size());
+    check(cudaMemcpy(actual.data(),device,actual.size()*sizeof(float),cudaMemcpyDeviceToHost));
+    for (size_t i=0;i<actual.size();++i)
+      if (!std::isfinite(actual[i]) || actual[i] != expected[i])
+        throw std::runtime_error(std::string(variantNames[variant])+" output or coverage mismatch at query "+std::to_string(i));
+    std::cout << "variant=" << variantNames[variant] << " phase=" << phase
+              << " queries_checked=" << actual.size() << " output_fnv64="
+              << std::hex << outputFingerprint(actual) << std::dec << " correctness=PASS\n";
+  }
+
+  void initializeResults(float *results, uint32_t count)
+  {
+    check(cudaMemset(results,0xff,size_t(count)*sizeof(float)));
+  }
+
+  void validateRootCases()
+  {
+    const Grid grid{3,2,32,100};
+    const auto dispatch = makeDispatch(grid);
+    DeviceBuffer<uint32_t> deviceDispatch(dispatch.size()), primitiveIDs(1);
+    DeviceBuffer<box2f> deviceBoxes(1);
+    DeviceBuffer<bvh2f::node_t> deviceNodes(1);
+    DeviceBuffer<float> output(dispatch.size());
+    check(cudaMemcpy(deviceDispatch.data,dispatch.data(),dispatch.size()*sizeof(uint32_t),cudaMemcpyHostToDevice));
+    const box2f box{cuBQL::vec2f{-120,-120},cuBQL::vec2f{-96,-96}};
+    check(cudaMemcpy(deviceBoxes.data,&box,sizeof(box),cudaMemcpyHostToDevice));
+    const uint32_t id = 0;
+    check(cudaMemcpy(primitiveIDs.data,&id,sizeof(id),cudaMemcpyHostToDevice));
+    bvh2f::node_t node{};
+    node.bounds = box;
+    node.admin.offset = 0;
+    node.admin.count = 1;
+    check(cudaMemcpy(deviceNodes.data,&node,sizeof(node),cudaMemcpyHostToDevice));
+    bvh2f single;
+    single.nodes = deviceNodes.data;
+    single.numNodes = 1;
+    single.primIDs = primitiveIDs.data;
+    single.numPrims = 1;
+    const std::vector<box2f> boxes{box};
+    std::vector<float> expected(dispatch.size(),0);
+    for (uint32_t query=0;query<expected.size();++query)
+      expected[query] = float(bruteQuery(query,boxes,grid).area);
+    for (int variant=0;variant<3;++variant) {
+      initializeResults(output.data,uint32_t(dispatch.size()));
+      launch(variant,single,single,deviceBoxes.data,deviceBoxes.data,grid,deviceDispatch.data,output.data);
+      check(cudaDeviceSynchronize());
+      validate(expected,output.data,variant,"single_root_leaf");
+      initializeResults(output.data,uint32_t(dispatch.size()));
+      launch(variant,bvh2f{},bvh2f{},nullptr,nullptr,grid,deviceDispatch.data,output.data);
+      check(cudaDeviceSynchronize());
+      validate(std::vector<float>(dispatch.size(),0),output.data,variant,"empty_tree");
+    }
+  }
+
+  float measure(int variant, bvh2f reference, bvh2f optimized, const box2f *boxes,
+                const box2f *reordered, Grid grid, const uint32_t *dispatch,
+                float *results, Event &start, Event &stop)
+  {
+    initializeResults(results,grid.x*grid.y);
+    check(cudaEventRecord(start.handle));
+    launch(variant,reference,optimized,boxes,reordered,grid,dispatch,results);
+    check(cudaEventRecord(stop.handle));
+    check(cudaEventSynchronize(stop.handle));
+    float milliseconds = 0;
+    check(cudaEventElapsedTime(&milliseconds,start.handle,stop.handle));
+    if (!std::isfinite(milliseconds) || milliseconds <= 0)
+      throw std::runtime_error("invalid CUDA event interval");
+    return milliseconds;
+  }
+
+  double median(std::vector<float> samples)
+  {
+    std::sort(samples.begin(),samples.end());
+    return 0.5*(double(samples[samples.size()/2-1])+samples[samples.size()/2]);
+  }
+
+  void run(const Options &options)
+  {
+    const Grid grid{options.queriesX,options.queriesY,options.stride,options.window};
+    const uint32_t count = grid.x*grid.y;
+    int deviceID = 0, runtimeVersion = 0, driverVersion = 0;
+    cudaDeviceProp device{};
+    check(cudaGetDevice(&deviceID));
+    check(cudaGetDeviceProperties(&device,deviceID));
+    check(cudaRuntimeGetVersion(&runtimeVersion));
+    check(cudaDriverGetVersion(&driverVersion));
+    if (device.warpSize != warpUnified::warpSize || device.maxThreadsPerBlock < referenceBlockSize
+        || uint64_t(count+1)/2 > uint64_t(device.maxGridSize[0]))
+      throw std::runtime_error("device does not support the required warp or launch geometry");
+    std::cout << std::fixed << std::setprecision(6) << "device=" << device.name
+              << " runtime=" << runtimeVersion << " driver=" << driverVersion << '\n';
+    validateRootCases();
+    const auto boxes = eda_layout::makeBoxes<float,2>(options.tilesX,options.tilesY,
+                                                    eda_layout::Pattern::Layout);
+    std::cout << "input=shared_morton_layout tiles=" << options.tilesX << 'x' << options.tilesY
+              << " boxes=" << boxes.size() << " input_bytes=" << boxes.size()*sizeof(box2f)
+              << " geometry_fnv64=" << std::hex << eda_layout::fingerprint(boxes) << std::dec << '\n';
+    std::cout << "queries=" << grid.x << 'x' << grid.y << " query_count=" << count
+              << " origin=" << grid.originX << ',' << grid.originY << " stride=" << grid.stride
+              << " window=" << grid.window << " windows_overlap=" << (grid.window > grid.stride ? "yes" : "no") << '\n';
+    std::cout << "reference=radixBuilder leaf=1 max_allowed_leaf=1 block_threads=" << referenceBlockSize
+              << " queries_per_thread=1 query_order=row_major box_lookup=primitive_ids\n"
+              << "before=radixBuilder leaf=16 block_threads=" << blockSize
+              << " queries_per_warp=1 query_order=morton box_lookup=leaf_order\n"
+              << "after=radixBuilder leaf=16 block_threads=" << blockSize
+              << " queries_per_warp=1 query_order=morton box_lookup=leaf_order\n";
+    const auto expected = expectedResults(options,boxes,grid);
+    const auto dispatch = makeDispatch(grid);
+    DeviceBuffer<box2f> deviceBoxes(boxes.size()), reordered(boxes.size());
+    DeviceBuffer<uint32_t> deviceDispatch(dispatch.size());
+    DeviceBuffer<float> reference(count), before(count), after(count);
+    float *outputs[] = {reference.data,before.data,after.data};
+    check(cudaMemcpy(deviceBoxes.data,boxes.data(),boxes.size()*sizeof(box2f),cudaMemcpyHostToDevice));
+    check(cudaMemcpy(deviceDispatch.data,dispatch.data(),dispatch.size()*sizeof(uint32_t),cudaMemcpyHostToDevice));
+    Tree referenceTree, tree;
+    auto referenceConfig = cuBQL::BuildConfig(1);
+    referenceConfig.maxAllowedLeafSize = 1;
+    cuBQL::cuda::radixBuilder(referenceTree.bvh,deviceBoxes.data,uint32_t(boxes.size()),
+                             referenceConfig,0,referenceTree.resource);
+    cuBQL::cuda::radixBuilder(tree.bvh,deviceBoxes.data,uint32_t(boxes.size()),
+                             cuBQL::BuildConfig(16),0,tree.resource);
+    warpUnified::gatherBoxesLeafOrderKernel<<<uint32_t((boxes.size()+blockSize-1)/blockSize),blockSize>>>(
+      deviceBoxes.data,tree.bvh.primIDs,int(boxes.size()),reordered.data);
+    check(cudaGetLastError());
+    check(cudaDeviceSynchronize());
+    std::cout << "variant=reference ";
+    auditTree(referenceTree.bvh,boxes,nullptr,referenceStackDepth);
+    std::cout << "variant=before_after ";
+    auditTree(tree.bvh,boxes,reordered.data,warpUnified::stackDepth);
+    for (int variant=0;variant<3;++variant) {
+      for (int warmup=0;warmup<(options.benchmark ? warmupCount : 1);++warmup) {
+        initializeResults(outputs[variant],count);
+        launch(variant,referenceTree.bvh,tree.bvh,deviceBoxes.data,reordered.data,
+               grid,deviceDispatch.data,outputs[variant]);
+        check(cudaDeviceSynchronize());
+      }
+      validate(expected,outputs[variant],variant,"before_timing");
+    }
+    if (!options.benchmark) {
+      std::cout << "status=PASS mode=correctness\n";
+      return;
+    }
+    Event start, stop;
+    constexpr int order[6][3] = {{0,1,2},{0,2,1},{1,0,2},{1,2,0},{2,0,1},{2,1,0}};
+    std::vector<float> samples[3];
+    for (int round=0;round<roundCount;++round)
+      for (int position=0;position<3;++position) {
+        const int variant = order[round%6][position];
+        const float milliseconds = measure(variant,referenceTree.bvh,tree.bvh,deviceBoxes.data,
+                                           reordered.data,grid,deviceDispatch.data,outputs[variant],start,stop);
+        samples[variant].push_back(milliseconds);
+        std::cout << "round=" << round+1 << " position=" << position+1
+                  << " variant=" << variantNames[variant] << " query_kernel_ms=" << milliseconds << '\n';
+      }
+    for (int variant=0;variant<3;++variant)
+      validate(expected,outputs[variant],variant,"after_timing");
+    const double referenceMedian = median(samples[0]);
+    const double beforeMedian = median(samples[1]), afterMedian = median(samples[2]);
+    std::cout << "median_reference_ms=" << referenceMedian << " median_before_ms=" << beforeMedian
+              << " median_after_ms=" << afterMedian << " reference_to_andrea=" << referenceMedian/beforeMedian
+              << " reference_to_lane0=" << referenceMedian/afterMedian << " andrea_to_lane0=" << beforeMedian/afterMedian
+              << " reference_to_andrea_reduction_pct=" << 100*(1-beforeMedian/referenceMedian)
+              << " reference_to_lane0_reduction_pct=" << 100*(1-afterMedian/referenceMedian)
+              << " andrea_to_lane0_reduction_pct=" << 100*(1-afterMedian/beforeMedian) << '\n';
+    std::cout << "status=PASS timing_scope=query_kernel_only balanced_triplets=" << roundCount
+              << " warmups_per_variant=" << warmupCount
+              << " reference=scalar_row_major_radix_leaf1_block128"
+                 " before=bounds_broadcast_all_lane_overlap after=lane0_overlap_broadcast\n";
+  }
+}
+
+int main(int argc, char **argv)
+{
+  try {
+    const auto options = parseOptions(argc,argv);
+    if (options.help)
+      printHelp();
+    else
+      run(options);
+    return 0;
+  } catch (const std::exception &error) {
+    std::cerr << "warp_unified failed: " << error.what() << '\n';
+    return 1;
+  }
+}
diff --git a/tests/warp_unified.cuh b/tests/warp_unified.cuh
new file mode 100644
index 0000000..23122ff
--- /dev/null
+++ b/tests/warp_unified.cuh
@@ -0,0 +1,191 @@
+// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+// SPDX-License-Identifier: Apache-2.0
+
+#pragma once
+
+// Adapted from Andrea Paris' WarpUnified traversal (2a0f4db) and
+// the lane-0 child-overlap change (96b9db49).
+
+#include "cuBQL/bvh.h"
+#include <cuda_runtime.h>
+
+namespace warpUnified {
+
+  using bvh2f = cuBQL::BinaryBVH<float,2>;
+  using box2f = cuBQL::box_t<float,2>;
+  using Admin = bvh2f::node_t::Admin;
+
+  constexpr int warpSize = 32;
+  constexpr int stackDepth = 32;
+  constexpr unsigned fullWarp = 0xffffffffu;
+
+  static_assert(sizeof(Admin) == sizeof(uint64_t),
+                "WarpUnified requires an 8-byte Admin");
+
+  inline __device__ Admin loadAdmin(const bvh2f::node_t *node)
+  {
+    Admin admin;
+    admin.offsetAndCountBits = __ldg(&node->admin.offsetAndCountBits);
+    return admin;
+  }
+
+  inline __device__ box2f loadBox(const box2f *box)
+  {
+    box2f result;
+    result.lower.x = __ldg(&box->lower.x);
+    result.lower.y = __ldg(&box->lower.y);
+    result.upper.x = __ldg(&box->upper.x);
+    result.upper.y = __ldg(&box->upper.y);
+    return result;
+  }
+
+  inline __device__ void broadcastBox(box2f &box)
+  {
+    box.lower.x = __shfl_sync(fullWarp,box.lower.x,0);
+    box.lower.y = __shfl_sync(fullWarp,box.lower.y,0);
+    box.upper.x = __shfl_sync(fullWarp,box.upper.x,0);
+    box.upper.y = __shfl_sync(fullWarp,box.upper.y,0);
+  }
+
+  inline __device__ bool overlaps(const box2f &query, const box2f &box)
+  {
+    return (query.lower.x <= box.upper.x)
+      & (query.lower.y <= box.upper.y)
+      & (box.lower.x <= query.upper.x)
+      & (box.lower.y <= query.upper.y);
+  }
+
+  inline __device__ float intersectionArea(const box2f &query,
+                                         const box2f &box)
+  {
+    const float x0 = fmaxf(query.lower.x,box.lower.x);
+    const float y0 = fmaxf(query.lower.y,box.lower.y);
+    const float x1 = fminf(query.upper.x,box.upper.x);
+    const float y1 = fminf(query.upper.y,box.upper.y);
+    const float dx = fmaxf(0.f,x1-x0);
+    const float dy = fmaxf(0.f,y1-y0);
+    return fmaf(dx,dy,0.f);
+  }
+
+  inline __device__ float accumulateLeafArea(const box2f &query,
+                                           const box2f *reordered,
+                                           uint32_t offset, uint32_t count)
+  {
+    const int lane = threadIdx.x & (warpSize-1);
+    float area = 0.f;
+    for (uint32_t i=lane;i<count;i+=warpSize)
+      area += intersectionArea(query,loadBox(reordered+offset+i));
+#pragma unroll
+    for (int offset=warpSize/2;offset>0;offset>>=1)
+      area += __shfl_down_sync(fullWarp,area,offset);
+    return area;
+  }
+
+  struct TraversalState {
+    int depth = 0;
+    Admin admin{};
+    bool active = true;
+  };
+
+  inline __device__ void pop(TraversalState &state, const Admin *stack)
+  {
+    if (state.depth > 0)
+      state.admin = stack[--state.depth];
+    else
+      state.active = false;
+  }
+
+  inline __device__ void broadcastState(TraversalState &state)
+  {
+    state.admin.offsetAndCountBits
+      = __shfl_sync(fullWarp,state.admin.offsetAndCountBits,0);
+    state.depth = __shfl_sync(fullWarp,state.depth,0);
+    state.active = __shfl_sync(fullWarp,int(state.active),0) != 0;
+  }
+
+  inline __device__ void descend(TraversalState &state, Admin *stack,
+                                const bvh2f &bvh, uint32_t pairBase,
+                                unsigned hits)
+  {
+    if (hits & 1u) {
+      if (hits & 2u) {
+        if (state.depth == stackDepth)
+          asm volatile("trap;");
+        stack[state.depth++] = loadAdmin(bvh.nodes+pairBase+1u);
+      }
+      state.admin = loadAdmin(bvh.nodes+pairBase);
+    } else if (hits & 2u) {
+      state.admin = loadAdmin(bvh.nodes+pairBase+1u);
+    } else {
+      state.admin.count = 0;
+      pop(state,stack);
+    }
+  }
+
+  // Called by a full warp with one query and a per-warp stack; lane 0 returns the area.
+  template<bool Lane0>
+  inline __device__ float queryArea(const bvh2f &bvh, box2f query,
+                                   const box2f *reordered, Admin *stack)
+  {
+    if (bvh.numPrims == 0)
+      return 0.f;
+
+    const int lane = threadIdx.x & (warpSize-1);
+    TraversalState state;
+    if (lane == 0)
+      state.admin = loadAdmin(bvh.nodes);
+    broadcastState(state);
+
+    float area = 0.f;
+    while (state.active) {
+      if (state.admin.count == 0) {
+        const uint32_t pairBase = uint32_t(state.admin.offset);
+        if constexpr (Lane0) {
+          if (lane == 0) {
+            const box2f bounds0 = loadBox(&bvh.nodes[pairBase].bounds);
+            const box2f bounds1 = loadBox(&bvh.nodes[pairBase+1u].bounds);
+            const unsigned hits = (overlaps(query,bounds0) ? 1u : 0u)
+              | (overlaps(query,bounds1) ? 2u : 0u);
+            descend(state,stack,bvh,pairBase,hits);
+          }
+        } else {
+          box2f bounds0;
+          box2f bounds1;
+          if (lane == 0) {
+            bounds0 = loadBox(&bvh.nodes[pairBase].bounds);
+            bounds1 = loadBox(&bvh.nodes[pairBase+1u].bounds);
+          }
+          broadcastBox(bounds0);
+          broadcastBox(bounds1);
+          const unsigned hits = (overlaps(query,bounds0) ? 1u : 0u)
+            | (overlaps(query,bounds1) ? 2u : 0u);
+          if (lane == 0)
+            descend(state,stack,bvh,pairBase,hits);
+        }
+        broadcastState(state);
+      }
+
+      if (state.active && state.admin.count != 0) {
+        const float leafArea
+          = accumulateLeafArea(query,reordered,uint32_t(state.admin.offset),
+                               uint32_t(state.admin.count));
+        if (lane == 0) {
+          area += leafArea;
+          pop(state,stack);
+        }
+        broadcastState(state);
+      }
+    }
+    return area;
+  }
+
+  __global__ void gatherBoxesLeafOrderKernel(const box2f *src,
+                                           const uint32_t *primIDs,
+                                           int count, box2f *dst)
+  {
+    const uint32_t i = blockIdx.x*blockDim.x+threadIdx.x;
+    if (i < uint32_t(count))
+      dst[i] = src[primIDs[i]];
+  }
+
+}