packed_half_tree_dpf.cuh

Go to the documentation of this file.

// SPDX-License-Identifier: Apache-2.0
 
#pragma once
#include <cuda_runtime.h>
#include <type_traits>
#include <cstddef>
#include <cassert>
#include <vector>
#ifdef _OPENMP
#include <omp.h>
#endif
#include <fss/hash.cuh>
#include <fss/util.cuh>
#include <fss/constraints.cuh>
#include <fss/compat.cuh>
#include <cuda/std/span>
#include <limits>
 
#include <fss/half_tree_dpf.cuh>
#include <stdexcept>
 
namespace fss {
 
namespace detail {
class SchemeAccess;
 
// Split Cws view shared by the host and device split paths. The packed
// construction carries no control bits, so layer l is just cws[l*stride] and the
// terminal entry holds the packed-output correction word at cws[kDepth*stride].
// operator[] recovers the Cw the shared algorithms consume, so their bodies stay
// identical on the Cw and split paths.
template <typename Cw>
struct PackedHalfTreeCwsView {
  const int4 *cws;
  size_t stride;
 
  __host__ __device__ inline Cw operator[](int l) const {
    return Cw{cws[static_cast<size_t>(l) * stride]};
  }
};
 
}  // namespace detail
 
template <int in_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> && CcrHashable<Hash>)
class HalfTreeGrottoDcf;

template <int in_bits, int out_bits, typename Hash, typename In = uint>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
class PackedHalfTreeDpf {
public:
  static constexpr int kLanes = 128 / out_bits;

  static constexpr int kDepth = in_bits - [] {
    int r = 0;
    int v = kLanes;
    while (v > 1) {
      v >>= 1;
      ++r;
    }
    return r;
  }();
  static_assert(kDepth >= 1, "in_bits must exceed log2(kLanes) so that the tree has at least one level");

  struct __align__(16) Cw {
    int4 s;  
  };
  static_assert(sizeof(Cw) == 16);

  static constexpr int kCwCount = kDepth + 1;
 
  using Input = In;  

  __host__ __device__ explicit PackedHalfTreeDpf(Hash hash, int4 instance_seed);

  __host__ __device__ void Gen(cuda::std::span<int4> cws,
      cuda::std::span<const int4, 2> s0s, In a, uint64_t beta, size_t stride = 1);

  __host__ __device__ int4 Eval(bool b, int4 s0, cuda::std::span<const int4> cws,
      In x, size_t stride = 1);

  static __host__ __device__ uint64_t Extract(int4 block, In x);

  void EvalAll(bool b, int4 s0, cuda::std::span<const int4> cws, int4 ys[],
      int par_depth = 0, size_t stride = 1);

  void EvalRange(bool b, In first, In last, int4 s0, cuda::std::span<const int4> cws,
      cuda::std::span<int4> ys, int par_depth = 0, size_t stride = 1);

  template <int frontier_depth = 16, int threads_per_block = 256>
  __host__ cudaError_t EvalAllGpu(bool b, int4 s0, cuda::std::span<const int4> cws,
      cuda::std::span<int4> ys, cudaStream_t stream = 0, size_t stride = 1) const
    requires(kDepth < sizeof(size_t) * 8);

  template <int frontier_depth = 16, int threads_per_block = 256>
  __host__ cudaError_t EvalAllGpuBatch(bool b, cuda::std::span<const int4> s0s,
      cuda::std::span<const int4> cws, cuda::std::span<int4> ys,
      cudaStream_t stream = 0, size_t stride = 0) const
    requires(kDepth < sizeof(size_t) * 8);
 
#if FSS_ENABLE_V1_COMPAT
  // Deprecated v1.3.0 aggregate-Key entries and public constants; see
  // fss/compat.cuh. The Key alias exposes the compat record under the legacy
  // nested name, while the method declarations spell the record type so a
  // deprecated name never enters a member signature. The constants restore the
  // legacy public names over the renamed private spellings; kTopDepth was a
  // constant of the old implementation and does not steer the current
  // scheduling. EvalAll() runs sequentially, as the current split-component
  // default does; the old automatic thread-count split is not restored.

  using Key [[deprecated("use the split-component Gen/Eval/EvalAll with the packed-output terminal in the cws array")]] =
      compat::detail::FcwKeyRecord<PackedHalfTreeDpf, Cw, kDepth>;

  [[deprecated("beta must fit in out_bits bits; express the bound from out_bits instead")]]
  static constexpr uint64_t kBetaMax = out_bits == 64 ? ~0ULL : (1ULL << out_bits) - 1;

  [[deprecated("the constant described the removed scalar-top implementation and does not steer the current scheduling")]]
  static constexpr int kTopDepth = 3;

  [[deprecated("use the split-component Gen(cuda::std::span<int4>, cuda::std::span<const int4, 2>, In, uint64_t, size_t) overload")]]
  __host__ __device__ void Gen(compat::detail::FcwKeyRecord<PackedHalfTreeDpf, Cw, kDepth> &k0,
      compat::detail::FcwKeyRecord<PackedHalfTreeDpf, Cw, kDepth> &k1,
      cuda::std::span<const int4, 2> s0s, In a, uint64_t beta);

  [[deprecated("use the split-component Eval(bool, int4, cuda::std::span<const int4>, In, size_t) overload")]]
  __host__ __device__ int4 Eval(bool b, const compat::detail::FcwKeyRecord<PackedHalfTreeDpf, Cw, kDepth> &key, In x);

  [[deprecated("use the split-component EvalAll(bool, int4, cuda::std::span<const int4>, int4[], int, size_t) overload")]]
  void EvalAll(bool b, const compat::detail::FcwKeyRecord<PackedHalfTreeDpf, Cw, kDepth> &key,
      cuda::std::span<int4> ys);
#endif  // FSS_ENABLE_V1_COMPAT
 
private:
  friend class detail::SchemeAccess;
 
  template <int frontier_depth, int threads_per_block>
  struct GpuLaunchShape;
 
  static __host__ __device__ void ExpandNode(Hash &hash, int4 instance_seed,
      int4 node, int4 cw, int4 &left, int4 &right);
  static __host__ __device__ int4 PathStep(Hash &hash, int4 instance_seed,
      int4 node, int4 cw, bool take_right);
  static __host__ __device__ int4 ConvertLeaf(Hash &hash, int4 instance_seed,
      int4 node, int4 output_cw);
 
  Hash hash_;
  int4 instance_seed_;
 
  // Split Cws view over the caller's array. The host methods and the
  // split GPU kernels share this one view; see detail::PackedHalfTreeCwsView.
  using CwsView = detail::PackedHalfTreeCwsView<Cw>;
 
  // Writes the strided Cw slots.
  struct CwsWriter {
    int4 *cws;
    size_t stride;
 
    __host__ __device__ void Write(int i, int4 s_cw) {
      cws[static_cast<size_t>(i) * stride] = s_cw;
    }
 
    __host__ __device__ void Finish(int4 terminal) {
      cws[static_cast<size_t>(kDepth) * stride] = terminal;
    }
  };
 
  // Shared Cw generation loop. The writer is passed by value and owns
  // only output access, so the existing algorithm and backend use stay unchanged.
  template <typename Writer>
  __host__ __device__ void Generate(Writer writer, const int4 *s0s, In a, uint64_t beta);
 
  // Shared point evaluation over Cw arrays and split Cws views.
  template <typename Cws>
  __host__ __device__ int4 EvalCws(bool b, int4 s0, Cws cws, In x);
 
  // Shared full-domain evaluation over Cw arrays and split Cws views.
  template <typename Cws>
  void EvalAllCws(bool b, int4 s0, Cws cws, int4 ys[], int par_depth);
 
  // Shared interval evaluation. Cws is a Cw pointer on the array form
  // and a CwsView on the split form.
  template <typename Cws>
  void EvalRangeCws(bool b, In first, In last, int4 s0, Cws cws,
      cuda::std::span<int4> ys, int par_depth);
 
  template <typename Cws>
  void EvalRangeTree(int4 node, Cws cws, int4 ys[], In first,
      In last, In output_first, In low, In high, int level);
 
  template <int other_bits, typename OtherHash, typename OtherIn>
    requires(ValidForInputBits<OtherIn, other_bits> && CcrHashable<OtherHash>)
  friend class HalfTreeGrottoDcf;
 
  // The private spellings below back the deprecated public aliases kBetaMax
  // and kTopDepth of the compatibility section; the values and every use are
  // unchanged by the rename.
  static constexpr uint64_t kBetaLimit = out_bits == 64 ? ~0ULL : (1ULL << out_bits) - 1;
  // Three scalar levels produce the eight parents consumed by Expand8().
  static constexpr int kBatchTopDepth = 3;
 
  static __host__ __device__ int4 PackBeta(In a, uint64_t beta);
 
  // Walks x's path to its level-kDepth node. visit(t, x_bit) sees the control bit of every
  // expanded node and the branch x takes from it.
  template <typename Cws, typename Visit>
  __host__ __device__ int4 Walk(bool b, int4 s0, Cws cws, In x, Visit &&visit);
 
  // lhs ^ (rhs & mask) with mask being 0 or -1 in every lane, so the per-node correction words
  // controlled by the seed LSB apply branchlessly.
  static __host__ __device__ int4 MaskXor(int4 lhs, int4 rhs, int mask);
 
  // Expands one node with a scalar hash call: left = H_S(node) ^ (t ? cw : 0), right = left ^ node.
  __host__ __device__ void ExpandNode(int4 node, int4 cw, int4 &left, int4 &right);
 
  // Level-order full-domain evaluation over a batched frontier. The frontier of
  // level i lives in bufs[i & 1]; the final expansion fuses the leaf conversion and writes
  // each ys block exactly once, in the same order the DFS walk produces. The level-order
  // node index equals the in-order leaf index, so the ys layout is unchanged.
  template <typename Cws>
  void EvalAllFrontier(bool b, int4 s0, Cws cws, int4 fcw, int4 ys[], int par_depth);
 
  // Expands the 8 frontier parents src[8k..8k+8) into their 16 children with one batch.
  // Parent p has children 2p and 2p + 1, so the children of a batch interleave in pairs.
  void Expand8(const int4 *src, int4 *dst, size_t k, int4 cw);
 
  // Last-level Expand8() with the leaf conversion fused in: the children are the level
  // kDepth nodes, so they are hashed straight from registers into two more batches
  // (even and odd children), fcw applies through the child LSB mask, and each final block
  // is stored exactly once.
  void ExpandConvert8(const int4 *src, int4 *ys, size_t k, int4 cw, int4 fcw);
 
  __host__ __device__ int4 ConvertLeaf(int4 node, int4 fcw);
};
 
// -----------------------------------------------------------------------------
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::PackedHalfTreeDpf(Hash hash, int4 instance_seed)
    : hash_(hash), instance_seed_(instance_seed) {}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline int4 PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::PackBeta(
    In a, uint64_t beta) {
  int lane = static_cast<int>(a) & (kLanes - 1);
  int4 packed = {0, 0, 0, 0};
  if constexpr (out_bits == 64) {
    packed = lane == 0 ? int4{static_cast<int>(beta), static_cast<int>(beta >> 32), 0, 0}
                       : int4{0, 0, static_cast<int>(beta), static_cast<int>(beta >> 32)};
  } else {
    constexpr int lanes_per_word = 32 / out_bits;
    int word = lane / lanes_per_word;
    int off = (lane % lanes_per_word) * out_bits;
    unsigned int v = static_cast<unsigned int>(beta) << off;
    if (word == 0) packed.x = static_cast<int>(v);
    else if (word == 1) packed.y = static_cast<int>(v);
    else if (word == 2) packed.z = static_cast<int>(v);
    else packed.w = static_cast<int>(v);
  }
  return packed;
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline uint64_t PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::Extract(
    int4 block, In x) {
  int lane = static_cast<int>(x) & (kLanes - 1);
  if constexpr (out_bits == 64) {
    uint64_t lo = lane == 0 ? static_cast<unsigned int>(block.x) : static_cast<unsigned int>(block.z);
    uint64_t hi = lane == 0 ? static_cast<unsigned int>(block.y) : static_cast<unsigned int>(block.w);
    return lo | (hi << 32);
  } else {
    constexpr int lanes_per_word = 32 / out_bits;
    int word = lane / lanes_per_word;
    int off = (lane % lanes_per_word) * out_bits;
    unsigned int w = word == 0 ? static_cast<unsigned int>(block.x)
        : word == 1 ? static_cast<unsigned int>(block.y)
        : word == 2 ? static_cast<unsigned int>(block.z)
        : static_cast<unsigned int>(block.w);
    return (w >> off) & static_cast<unsigned int>((1ULL << out_bits) - 1);
  }
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::Gen(
    cuda::std::span<int4> cws, cuda::std::span<const int4, 2> s0s, In a, uint64_t beta,
    size_t stride) {
  util::RequireArg(stride > 0, "PackedHalfTreeDpf::Gen: stride must be positive");
  constexpr size_t kMaxElements = SIZE_MAX / sizeof(int4);
  util::RequireArg(stride <= (kMaxElements - 1) / static_cast<size_t>(kDepth),
      "PackedHalfTreeDpf::Gen: stride exceeds the correction array limit");
  const size_t required = static_cast<size_t>(kDepth) * stride + 1;
  util::RequireArg(cws.size() >= required,
      "PackedHalfTreeDpf::Gen: cws must hold kDepth*stride+1 elements");
  Generate(CwsWriter{cws.data(), stride}, s0s.data(), a, beta);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <typename Writer>
__host__ __device__ inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::Generate(
    Writer writer, const int4 *s0s, In a, uint64_t beta) {
  util::RequireArg(beta <= kBetaLimit, "PackedHalfTreeDpf::Gen: beta exceeds kBetaMax");
 
  int4 packed_beta = PackBeta(a, beta);
 
  // node1 starts with the control bit set, so the path difference delta keeps LSB 1 and
  // exactly one party has t=1 on the path at every level.
  int4 node0 = util::SetLsb(s0s[0], false);
  int4 node1 = util::SetLsb(s0s[1], true);
  int4 delta = util::Xor(node0, node1);  // LSB = 0^1 = 1
 
  for (int i = 0; i < kDepth; ++i) {
    int4 h0 = hash_.Hash(util::Xor(instance_seed_, node0));
    int4 h1 = hash_.Hash(util::Xor(instance_seed_, node1));
 
    bool a_bit = (a >> (in_bits - 1 - i)) & 1;
 
    // CW = h0 ^ h1 ^ (!a_bit ? delta : 0)
    // When a_bit=0 (go left): non-alpha is right, CW = h0^h1^delta makes right0=right1
    // When a_bit=1 (go right): non-alpha is left, CW = h0^h1 makes left0=left1
    int4 cw = util::Xor(h0, h1);
    if (!a_bit) cw = util::Xor(cw, delta);
 
    writer.Write(i, cw);
 
    bool t0 = util::GetLsb(node0);
    bool t1 = util::GetLsb(node1);
 
    // node_b = h_b ^ (a_bit ? node_b : 0) ^ (t_b ? cw : 0)
    int4 zero4 = {0, 0, 0, 0};
    int4 ab_mask0 = a_bit ? node0 : zero4;
    int4 ab_mask1 = a_bit ? node1 : zero4;
    int4 t0_mask = t0 ? cw : zero4;
    int4 t1_mask = t1 ? cw : zero4;
 
    node0 = util::Xor(util::Xor(h0, ab_mask0), t0_mask);
    node1 = util::Xor(util::Xor(h1, ab_mask1), t1_mask);
 
    // delta = node0 ^ node1 for the next level; stays the root seed XOR
    delta = util::Xor(node0, node1);
  }
 
  // The leaf hash masks the constant path difference, see the file-level details.
  int4 leaf0 = hash_.Hash(util::Xor(instance_seed_, node0));
  int4 leaf1 = hash_.Hash(util::Xor(instance_seed_, node1));
  int4 fcw = util::Xor(util::Xor(packed_beta, leaf0), leaf1);
  writer.Finish(fcw);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline int4 PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::Eval(
    bool b, int4 s0, cuda::std::span<const int4> cws, In x, size_t stride) {
  util::RequireArg(stride > 0, "PackedHalfTreeDpf::Eval: stride must be positive");
  constexpr size_t kMaxElements = SIZE_MAX / sizeof(int4);
  util::RequireArg(stride <= (kMaxElements - 1) / static_cast<size_t>(kDepth),
      "PackedHalfTreeDpf::Eval: stride exceeds the correction array limit");
  const size_t required = static_cast<size_t>(kDepth) * stride + 1;
  util::RequireArg(cws.size() >= required,
      "PackedHalfTreeDpf::Eval: cws must hold kDepth*stride+1 elements");
  return EvalCws(b, s0, CwsView{cws.data(), stride}, x);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <typename Cws>
__host__ __device__ inline int4 PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalCws(
    bool b, int4 s0, Cws cws, In x) {
  return ConvertLeaf(Walk(b, s0, cws, x, [](bool, bool) {}), cws[kDepth].s);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <typename Cws, typename Visit>
__host__ __device__ inline int4 PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::Walk(
    bool b, int4 s0, Cws cws, In x, Visit &&visit) {
  int4 node = util::SetLsb(s0, b);
 
  for (int i = 0; i < kDepth; ++i) {
    bool x_bit = (x >> (in_bits - 1 - i)) & 1;
    bool t = util::GetLsb(node);
    visit(t, x_bit);
 
    node = PathStep(hash_, instance_seed_, node, cws[i].s, x_bit);
  }
 
  return node;
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalAll(
    bool b, int4 s0, cuda::std::span<const int4> cws, int4 ys[], int par_depth,
    size_t stride) {
  if (stride == 0) throw std::invalid_argument("PackedHalfTreeDpf::EvalAll: stride must be positive");
  constexpr size_t kMaxElements = SIZE_MAX / sizeof(int4);
  if (stride > (kMaxElements - 1) / static_cast<size_t>(kDepth))
    throw std::invalid_argument("PackedHalfTreeDpf::EvalAll: stride exceeds the correction array limit");
  const size_t required = static_cast<size_t>(kDepth) * stride + 1;
  if (cws.size() < required)
    throw std::invalid_argument("PackedHalfTreeDpf::EvalAll: cws must hold kDepth*stride+1 elements");
  EvalAllCws(b, s0, CwsView{cws.data(), stride}, ys, par_depth);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <typename Cws>
inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalAllCws(
    bool b, int4 s0, Cws cws, int4 ys[], int par_depth) {
  assert(par_depth >= 0);
  static_assert(detail::FitsFullDomain<int4, kDepth>(),
      "packed output count and byte size must fit in size_t");
  int4 fcw = cws[kDepth].s;
  if constexpr (kDepth > kBatchTopDepth) {
    EvalAllFrontier(b, s0, cws, fcw, ys, par_depth);
  } else {
    // Domains too shallow to fill an 8-wide frontier evaluate one point per
    // block through the tested point walk.
    constexpr size_t blocks = size_t{1} << kDepth;
    fss::detail::ParallelFor(par_depth, blocks, [&](size_t j) {
      ys[j] = EvalCws(b, s0, cws, static_cast<In>(j * kLanes));
    });
  }
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalRange(bool b, In first,
    In last, int4 s0, cuda::std::span<const int4> cws, cuda::std::span<int4> ys,
    int par_depth, size_t stride) {
  if (stride == 0)
    throw std::invalid_argument("PackedHalfTreeDpf::EvalRange: stride must be positive");
  constexpr size_t kMaxElements = SIZE_MAX / sizeof(int4);
  if (stride > (kMaxElements - 1) / static_cast<size_t>(kDepth))
    throw std::invalid_argument("PackedHalfTreeDpf::EvalRange: stride exceeds the correction array limit");
  const size_t required = static_cast<size_t>(kDepth) * stride + 1;
  if (cws.size() < required || cws.data() == nullptr)
    throw std::invalid_argument("PackedHalfTreeDpf::EvalRange: cws must hold kDepth*stride+1 elements");
  if (ys.data() == nullptr)
    throw std::invalid_argument("PackedHalfTreeDpf::EvalRange: ys must not be null");
  EvalRangeCws(b, first, last, s0, CwsView{cws.data(), stride}, ys, par_depth);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <typename Cws>
inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalRangeCws(
    bool b, In first, In last, int4 s0, Cws cws,
    cuda::std::span<int4> ys, int par_depth) {
  assert(par_depth >= 0);
  In domain_last = static_cast<In>(~In{0});
  if constexpr (in_bits < sizeof(In) * 8) {
    domain_last = static_cast<In>((In{1} << in_bits) - 1);
  }
  if (first > last)
    throw std::invalid_argument("PackedHalfTreeDpf::EvalRange: first must not exceed last");
  if (last > domain_last)
    throw std::invalid_argument("PackedHalfTreeDpf::EvalRange: last must stay inside the input domain");
  const __uint128_t distance = static_cast<__uint128_t>(last) - first;
  const __uint128_t block_distance = distance / kLanes;
  // The caller controls the interval, so this is an argument check, not an
  // invariant: with a 128-bit input domain the block distance can exceed
  // size_t, and narrowing it first would shrink the output-size check.
  if (block_distance >= std::numeric_limits<size_t>::max() / sizeof(int4))
    throw std::invalid_argument("PackedHalfTreeDpf::EvalRange: the interval is too large for size_t");
  const size_t blocks = static_cast<size_t>(block_distance) + 1;
  if (ys.size() < blocks || ys.data() == nullptr)
    throw std::invalid_argument("PackedHalfTreeDpf::EvalRange: ys must hold the covered blocks");
 
  const size_t partitions = detail::ResolveParPartitions(par_depth, blocks);
  detail::ParallelFor(par_depth, partitions, [&](size_t partition) {
    const size_t begin = partition * (blocks / partitions) +
        (partition < blocks % partitions ? partition : blocks % partitions);
    const size_t end = begin + blocks / partitions + (partition < blocks % partitions);
    for (size_t index = begin; index < end; ++index) ys[index] = {0, 0, 0, 0};
    const __uint128_t begin_offset = static_cast<__uint128_t>(begin) * kLanes;
    const __uint128_t end_offset = end == blocks ? distance :
        static_cast<__uint128_t>(end) * kLanes - 1;
    EvalRangeTree(util::SetLsb(s0, b), cws, ys.data(),
        static_cast<In>(first + begin_offset), static_cast<In>(first + end_offset),
        first, In{0}, domain_last, 0);
  });
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <typename Cws>
inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalRangeTree(
    int4 node, Cws cws, int4 ys[], In first, In last,
    In output_first, In low, In high, int level) {
  if (level == kDepth) {
    const int4 block = ConvertLeaf(node, cws[kDepth].s);
    const In begin = first > low ? first : low;
    const In end = last < high ? last : high;
    const size_t lanes = static_cast<size_t>(end - begin) + 1;
    for (size_t lane = 0; lane < lanes; ++lane) {
      const In x = static_cast<In>(begin + lane);
      const In offset = static_cast<In>(x - output_first);
      const size_t index = static_cast<size_t>(offset / kLanes);
      ys[index] = util::Xor(ys[index], PackBeta(offset, Extract(block, x)));
    }
    return;
  }
  const In middle = static_cast<In>(low + ((high - low) >> 1));
  int4 left, right;
  ExpandNode(node, cws[level].s, left, right);
  if (first <= middle) {
    EvalRangeTree(left, cws, ys, first, last, output_first, low, middle, level + 1);
  }
  if (last > middle) {
    EvalRangeTree(right, cws, ys, first, last, output_first,
        static_cast<In>(middle + 1), high, level + 1);
  }
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline int4 PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::MaskXor(
    int4 lhs, int4 rhs, int mask) {
  return {lhs.x ^ (rhs.x & mask), lhs.y ^ (rhs.y & mask), lhs.z ^ (rhs.z & mask), lhs.w ^ (rhs.w & mask)};
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::ExpandNode(
    int4 node, int4 cw, int4 &left, int4 &right) {
  ExpandNode(hash_, instance_seed_, node, cw, left, right);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <typename Cws>
inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalAllFrontier(
    bool b, int4 s0, Cws cws, int4 fcw, int4 ys[], int par_depth) {
  constexpr size_t num_blocks = detail::FullDomainSize<int4, kDepth>();
 
  // Frontier ping-pong scratch, grown on demand and reused across calls of this thread.
  // Two halves of num_blocks / 2 each: no frontier before the last level exceeds that size.
  // thread_local isolates concurrent calls: each thread pays the first-call allocation and
  // retains the peak 16 * 2 ** kDepth bytes, about one packed output buffer, until it exits.
  static thread_local std::vector<int4> scratch;
  if (scratch.size() < num_blocks) scratch.resize(num_blocks);
  int4 *bufs[2] = {scratch.data(), scratch.data() + num_blocks / 2};
 
  // Scalar top: the root becomes the 2^kBatchTopDepth level-kBatchTopDepth frontier in bufs[kBatchTopDepth & 1].
  bufs[0][0] = util::SetLsb(s0, b);
  for (int i = 0; i < kBatchTopDepth; ++i) {
    size_t count = size_t{1} << i;
    for (size_t j = 0; j < count; ++j) {
      ExpandNode(bufs[i & 1][j], cws[i].s, bufs[(i + 1) & 1][2 * j], bufs[(i + 1) & 1][2 * j + 1]);
    }
  }
 
  // Each level bounds its independent partitions by 2^par_depth. The
  // sequential branch avoids an inactive OpenMP region as well as workers.
  int depth = fss::detail::ResolveParDepth(par_depth, kDepth);
  if (depth == 0) {
    for (int i = kBatchTopDepth; i < kDepth; ++i) {
      const int4 *src = bufs[i & 1];
      int4 *dst = i + 1 == kDepth ? ys : bufs[(i + 1) & 1];
      size_t batches = (size_t{1} << i) / 8;
      for (size_t k = 0; k < batches; ++k) {
        if (i + 1 == kDepth) ExpandConvert8(src, dst, k, cws[i].s, fcw);
        else Expand8(src, dst, k, cws[i].s);
      }
    }
    return;
  }
  int threads = fss::detail::ResolveParThreads(depth, num_blocks / 16);
#pragma omp parallel num_threads(threads)
  {
    for (int i = kBatchTopDepth; i < kDepth; ++i) {
      const int4 *src = bufs[i & 1];
      int4 *dst = i + 1 == kDepth ? ys : bufs[(i + 1) & 1];
      size_t batches = (size_t{1} << i) / 8;
      size_t partitions = fss::detail::ResolveParPartitions(depth, batches);
#pragma omp for schedule(static)
      for (size_t partition = 0; partition < partitions; ++partition) {
        size_t begin = partition * (batches / partitions) +
            (partition < batches % partitions ? partition : batches % partitions);
        size_t end = begin + batches / partitions + (partition < batches % partitions);
        for (size_t k = begin; k < end; ++k) {
          if (i + 1 == kDepth) ExpandConvert8(src, dst, k, cws[i].s, fcw);
          else Expand8(src, dst, k, cws[i].s);
        }
      }
    }
  }
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::Expand8(
    const int4 *src, int4 *dst, size_t k, int4 cw) {
  int4 parents[8];
  cuda::std::array<int4, 8> seeds;
  for (int j = 0; j < 8; ++j) {
    parents[j] = src[8 * k + j];
    seeds[j] = util::Xor(instance_seed_, parents[j]);
  }
  auto hs = detail::HashBatch(hash_, seeds);
  for (int j = 0; j < 8; ++j) {
    int4 left = MaskXor(hs[j], cw, -static_cast<int>(util::GetLsb(parents[j])));
    dst[16 * k + 2 * j] = left;
    dst[16 * k + 2 * j + 1] = util::Xor(left, parents[j]);
  }
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::ExpandConvert8(
    const int4 *src, int4 *ys, size_t k, int4 cw, int4 fcw) {
  int4 parents[8];
  cuda::std::array<int4, 8> seeds;
  for (int j = 0; j < 8; ++j) {
    parents[j] = src[8 * k + j];
    seeds[j] = util::Xor(instance_seed_, parents[j]);
  }
  auto hs = detail::HashBatch(hash_, seeds);
 
  int4 children[16];
  cuda::std::array<int4, 8> conv_seeds[2];
  for (int j = 0; j < 8; ++j) {
    int4 left = MaskXor(hs[j], cw, -static_cast<int>(util::GetLsb(parents[j])));
    children[2 * j] = left;
    children[2 * j + 1] = util::Xor(left, parents[j]);
    conv_seeds[0][j] = util::Xor(instance_seed_, children[2 * j]);
    conv_seeds[1][j] = util::Xor(instance_seed_, children[2 * j + 1]);
  }
  auto leaves_even = detail::HashBatch(hash_, conv_seeds[0]);
  auto leaves_odd = detail::HashBatch(hash_, conv_seeds[1]);
  for (int j = 0; j < 8; ++j) {
    ys[16 * k + 2 * j] = MaskXor(leaves_even[j], fcw, -static_cast<int>(util::GetLsb(children[2 * j])));
    ys[16 * k + 2 * j + 1] =
        MaskXor(leaves_odd[j], fcw, -static_cast<int>(util::GetLsb(children[2 * j + 1])));
  }
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline int4 PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::ConvertLeaf(
    int4 node, int4 fcw) {
  return ConvertLeaf(hash_, instance_seed_, node, fcw);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::ExpandNode(
    Hash &hash, int4 instance_seed, int4 node, int4 cw, int4 &left, int4 &right) {
  int4 hashed = hash.Hash(util::Xor(instance_seed, node));
  left = MaskXor(hashed, cw, -static_cast<int>(util::GetLsb(node)));
  right = util::Xor(left, node);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline int4 PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::PathStep(
    Hash &hash, int4 instance_seed, int4 node, int4 cw, bool take_right) {
  const int4 zero = {0, 0, 0, 0};
  int4 hashed = hash.Hash(util::Xor(instance_seed, node));
  return util::Xor(util::Xor(hashed, take_right ? node : zero),
      util::GetLsb(node) ? cw : zero);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline int4 PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::ConvertLeaf(
    Hash &hash, int4 instance_seed, int4 node, int4 output_cw) {
  int4 leaf = hash.Hash(util::Xor(instance_seed, node));
  if (util::GetLsb(node)) leaf = util::Xor(leaf, output_cw);
  return leaf;
}
 
namespace gpu::detail {
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
struct HalfTreeGpuTraits<PackedHalfTreeDpf<in_bits, out_bits, Hash, In>> {
  using Backend = Hash;
  static constexpr int kExpansionDepth = PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::kDepth;
  static constexpr int kOutputBits = PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::kDepth;
  static constexpr bool kPacked = true;
};
 
#if defined(__CUDACC__)
 
// Split-layout counterparts of the shared HalfTree kernels. The Cw
// buffer arrives type-erased and the boundary view reconstructs the per-layer
// Cw entries, so the shared traversal stays identical on both paths.
template <typename Scheme, int frontier_depth, int block_root_depth>
__global__ void PackedHalfTreeSplitEvalAllKernel(bool b, int4 s0, const void *cws_buffer,
    size_t stride, int4 *ys, int4 instance_seed, typename HalfTreeGpuTraits<Scheme>::Backend hash) {
  using Cws = fss::detail::PackedHalfTreeCwsView<typename Scheme::Cw>;
  const auto *cws = static_cast<const int4 *>(cws_buffer);
  const int64_t thread_index = static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
  if (thread_index >= (1LL << frontier_depth)) return;
  __shared__ int4 shared_frontier[1 << (frontier_depth - block_root_depth)];
  HalfTreeEvalAllBlock<Scheme, frontier_depth, block_root_depth>(b, s0,
      Cws{cws, stride}, ys, instance_seed, hash, blockIdx.x, thread_index,
      shared_frontier);
}
 
template <typename Scheme, int frontier_depth, int block_root_depth>
__global__ void PackedHalfTreeSplitEvalAllBatchKernel(bool b, const int4 *s0s,
    const void *cws_buffer, size_t stride, int4 *ys, int4 instance_seed,
    typename HalfTreeGpuTraits<Scheme>::Backend hash) {
  using Cws = fss::detail::PackedHalfTreeCwsView<typename Scheme::Cw>;
  const auto *cws = static_cast<const int4 *>(cws_buffer);
  constexpr int kBlocksPerKey = 1 << block_root_depth;
  constexpr size_t kOutputsPerKey = size_t{1} << HalfTreeGpuTraits<Scheme>::kOutputBits;
  const int key_index = blockIdx.x / kBlocksPerKey;
  const int subtree_block_index = blockIdx.x - key_index * kBlocksPerKey;
  const int64_t thread_index = static_cast<int64_t>(subtree_block_index) * blockDim.x + threadIdx.x;
  if (thread_index >= (1LL << frontier_depth)) return;
  int4 *key_outputs = ys + static_cast<size_t>(key_index) * kOutputsPerKey;
  __shared__ int4 shared_frontier[1 << (frontier_depth - block_root_depth)];
  HalfTreeEvalAllBlock<Scheme, frontier_depth, block_root_depth>(b, s0s[key_index],
      Cws{cws + key_index, stride}, key_outputs, instance_seed, hash,
      subtree_block_index, thread_index, shared_frontier);
}
 
// Key-major batch counterpart: key k's components occupy its whole block of
// records with a unit stride.
template <typename Scheme, int frontier_depth, int block_root_depth>
__global__ void PackedHalfTreeSplitEvalAllBatchDeviceKeyKernel(bool b, const int4 *s0s,
    const void *cws_buffer, size_t key_records, size_t key_count, int4 *ys, int4 instance_seed,
    typename HalfTreeGpuTraits<Scheme>::Backend hash) {
  using Cws = fss::detail::PackedHalfTreeCwsView<typename Scheme::Cw>;
  const auto *cws = static_cast<const int4 *>(cws_buffer);
  constexpr int kBlocksPerKey = 1 << block_root_depth;
  constexpr size_t kOutputsPerKey = size_t{1} << HalfTreeGpuTraits<Scheme>::kOutputBits;
  const int key_index = blockIdx.x / kBlocksPerKey;
  const int subtree_block_index = blockIdx.x - key_index * kBlocksPerKey;
  const int64_t thread_index = static_cast<int64_t>(subtree_block_index) * blockDim.x + threadIdx.x;
  if (thread_index >= (1LL << frontier_depth)) return;
  const int4 *key = cws + static_cast<size_t>(key_index) * key_records;
  int4 *key_outputs = ys + static_cast<size_t>(key_index) * kOutputsPerKey;
  __shared__ int4 shared_frontier[1 << (frontier_depth - block_root_depth)];
  HalfTreeEvalAllBlock<Scheme, frontier_depth, block_root_depth>(b, s0s[key_index],
      Cws{key, 1}, key_outputs, instance_seed, hash,
      subtree_block_index, thread_index, shared_frontier);
}
 
// Host launcher for the key-major batch kernel; the hash backend and instance
// seed arrive by value from a host mirror of the device scheme.
template <int in_bits, int out_bits, typename Hash, typename In, int frontier_depth = 16,
    int threads_per_block = 256>
__host__ inline cudaError_t PackedHalfTreeDpfEvalAllGpuBatchDeviceKey(bool b, const int4 *s0s,
    size_t key_count, const int4 *cws, size_t cws_count, size_t key_records, int4 *ys,
    size_t ys_count, int4 instance_seed, const Hash &hash, cudaStream_t stream) {
  using Scheme = PackedHalfTreeDpf<in_bits, out_bits, Hash, In>;
  constexpr int kDepth = Scheme::kDepth;
  constexpr size_t kOutputsPerKey = size_t{1} << kDepth;
  constexpr size_t kMaxElements = std::numeric_limits<size_t>::max() / sizeof(int4);
  constexpr size_t kMaxAddressableSize = std::numeric_limits<size_t>::max();
  constexpr int kFrontierDepth = frontier_depth < kDepth ? frontier_depth : kDepth;
  static_assert(threads_per_block >= 1 && threads_per_block <= 1024 &&
      (threads_per_block & (threads_per_block - 1)) == 0,
      "gpu block size must be a power of two in [1,1024]");
  static_assert(kDepth - kFrontierDepth <= 8,
      "gpu frontier leaves more than eight tail levels; specify a larger frontier_depth");
  constexpr int kThreadsPerBlockLog2 = [] {
    int requested_log2 = 0;
    for (int threads = threads_per_block; threads > 1; threads >>= 1) ++requested_log2;
    return requested_log2 < kFrontierDepth ? requested_log2 : kFrontierDepth;
  }();
  constexpr int kBlockRootDepth = kFrontierDepth - kThreadsPerBlockLog2;
  constexpr int kBlocksPerKey = 1 << kBlockRootDepth;
  if (key_count == 0) return cudaSuccess;
  if (key_count > kMaxElements ||
      key_count > kMaxAddressableSize / kOutputsPerKey / sizeof(int4) ||
      key_count > kMaxAddressableSize / key_records / sizeof(int4) ||
      cws_count < key_count * key_records || ys_count < key_count * kOutputsPerKey ||
      s0s == nullptr || cws == nullptr || ys == nullptr) return cudaErrorInvalidValue;
  if (key_count > static_cast<size_t>(std::numeric_limits<int>::max()) / kBlocksPerKey)
    return cudaErrorInvalidConfiguration;
  PackedHalfTreeSplitEvalAllBatchDeviceKeyKernel<Scheme, kFrontierDepth, kBlockRootDepth>
      <<<key_count * kBlocksPerKey, 1 << kThreadsPerBlockLog2, 0, stream>>>(
          b, s0s, cws, key_records, key_count, ys, instance_seed, hash);
  return cudaGetLastError();
}
 
#endif  // defined(__CUDACC__)
 
#if defined(__CUDACC__)
 
// Point-evaluate packed HalfTree keys, one thread per key: ys[tid] is x's output
// block, or with cmp the 1[alpha <= x] share of HalfTreeGrottoDcf::Eval(). cws
// holds dense level-major correction words, the layout the split form of Gen()
// writes: layer i of key k is cws[i * nkeys + k] and the final correction word
// is the tail row, cws[kDepth * nkeys + k].
template <int in_bits, int unroll, bool b, bool cmp, typename Scheme, typename In, typename Query, typename Y>
__global__ void PackedHalfTreeDpfEvalPointKernel(Y *ys, const int4 *seeds, const int4 *cws, const Query *xs,
    int nkeys, typename HalfTreeGpuTraits<Scheme>::Backend hash, int4 instance_seed) {
  int tid = blockIdx.x * blockDim.x + threadIdx.x;
  if (tid >= nkeys) return;
  int4 node = fss::util::SetLsb(seeds[tid], b);
  In x = static_cast<In>(xs[tid]);
  fss::detail::GrottoTurns turns;
 
#pragma unroll unroll
  for (int i = 0; i < Scheme::kDepth; ++i) {
    bool x_bit = GetBit(x, in_bits - 1 - i);
    bool t = fss::util::GetLsb(node);
    if constexpr (cmp) turns(t, x_bit);
    int4 h = hash.Hash(fss::util::Xor(instance_seed, node));
    int4 cw = cws[(size_t)i * nkeys + tid];
    node = fss::util::Xor(fss::util::Xor(h, x_bit ? node : Zero4()), t ? cw : Zero4());
  }
 
  int4 leaf = hash.Hash(fss::util::Xor(instance_seed, node));
  if (fss::util::GetLsb(node)) leaf = fss::util::Xor(leaf, cws[(size_t)Scheme::kDepth * nkeys + tid]);
  if constexpr (cmp) {
    ys[tid] = turns.share ^ fss::detail::SchemeAccess::HalfTreeGrottoLeafParity<Scheme>(leaf, x);
  } else {
    ys[tid] = leaf;
  }
}
 
template <int unroll, bool cmp, int in_bits, typename Scheme, typename In, typename Query, typename Y>
void PackedHalfTreeDpfEvalPointLaunch(bool b, const int4 *seeds, const int4 *cws, const Query *xs, Y *ys,
    int nkeys, typename HalfTreeGpuTraits<Scheme>::Backend hash, int4 instance_seed, cudaStream_t stream) {
  if (nkeys == 0) return;
  constexpr int kBs = 256;
  int blocks = (nkeys + kBs - 1) / kBs;
  if (b) {
    PackedHalfTreeDpfEvalPointKernel<in_bits, unroll, true, cmp, Scheme, In, Query, Y>
        <<<blocks, kBs, 0, stream>>>(ys, seeds, cws, xs, nkeys, hash, instance_seed);
  } else {
    PackedHalfTreeDpfEvalPointKernel<in_bits, unroll, false, cmp, Scheme, In, Query, Y>
        <<<blocks, kBs, 0, stream>>>(ys, seeds, cws, xs, nkeys, hash, instance_seed);
  }
}

template <int unroll = 4, int in_bits, int out_bits, typename Hash, typename In, typename Query = In>
  requires QueryFor<Query, in_bits>
void PackedHalfTreeDpfEvalPointGpu(bool b, cuda::std::span<const int4> seeds, cuda::std::span<const int4> cws,
    cuda::std::span<const Query> xs, int4 *ys, int nkeys, const PackedHalfTreeDpf<in_bits, out_bits, Hash, In> &dpf,
    cudaStream_t stream = 0) {
  if (nkeys < 0)
    throw std::invalid_argument("PackedHalfTreeDpfEvalPointGpu: nkeys must be nonnegative");
  if (nkeys == 0) return;
  const size_t count = static_cast<size_t>(nkeys);
  constexpr size_t kLevels = static_cast<size_t>(PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::kDepth) + 1;
  if (count > SIZE_MAX / kLevels || cws.size() < count * kLevels)
    throw std::invalid_argument("PackedHalfTreeDpfEvalPointGpu: cws must hold (kDepth+1)*nkeys elements");
  if (seeds.size() < count || xs.size() < count || ys == nullptr)
    throw std::invalid_argument("PackedHalfTreeDpfEvalPointGpu: seeds, xs, and ys must cover nkeys entries");
  using Scheme = PackedHalfTreeDpf<in_bits, out_bits, Hash, In>;
  PackedHalfTreeDpfEvalPointLaunch<unroll, false, in_bits, Scheme, In, Query, int4>(b, seeds.data(), cws.data(),
      xs.data(), ys, nkeys, fss::detail::SchemeAccess::GetHash(dpf),
      fss::detail::SchemeAccess::GetInstanceSeed(dpf), stream);
}
 
#endif  // defined(__CUDACC__)
 
}  // namespace gpu::detail
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <int frontier_depth, int threads_per_block>
struct PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::GpuLaunchShape {
  static_assert(frontier_depth >= 0, "gpu frontier depth must be nonnegative");
  static_assert(1 <= threads_per_block && threads_per_block <= 1024 &&
      (threads_per_block & (threads_per_block - 1)) == 0,
      "gpu block size must be a power of two in [1,1024]");
  static_assert(kDepth >= 0, "gpu tree must have a nonnegative expansion depth");
  static constexpr int kFrontierDepth = frontier_depth < kDepth ? frontier_depth : kDepth;
  static_assert(kDepth - kFrontierDepth <= 8,
      "gpu frontier leaves more than eight tail levels; specify a larger frontier_depth");
  static constexpr int kThreadsPerBlockLog2 = [] {
    int depth = 0;
    for (int count = threads_per_block; count > 1; count >>= 1) ++depth;
    return depth < kFrontierDepth ? depth : kFrontierDepth;
  }();
  static constexpr int kThreadsPerBlock = 1 << kThreadsPerBlockLog2;
  static constexpr int kBlockRootDepth = kFrontierDepth - kThreadsPerBlockLog2;
  static_assert(kBlockRootDepth <= 30, "gpu block index exceeds the signed-int range");
  static constexpr size_t kBlocksPerKey = size_t{1} << (kBlockRootDepth <= 30 ? kBlockRootDepth : 0);
};
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <int frontier_depth, int threads_per_block>
__host__ inline cudaError_t PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalAllGpu(
    bool b, int4 s0, cuda::std::span<const int4> cws, cuda::std::span<int4> ys,
    cudaStream_t stream, size_t stride) const
  requires(kDepth < sizeof(size_t) * 8) {
#if defined(__CUDACC__)
  using LaunchShape = GpuLaunchShape<frontier_depth, threads_per_block>;
  using Scheme = PackedHalfTreeDpf<in_bits, out_bits, Hash, In>;
  constexpr size_t kOutputsPerKey = size_t{1} << kDepth;
  constexpr size_t kMaxElements = std::numeric_limits<size_t>::max() / sizeof(int4);
  if (stride == 0 || stride > (kMaxElements - 1) / static_cast<size_t>(kDepth))
    return cudaErrorInvalidValue;
  const size_t required = static_cast<size_t>(kDepth) * stride + 1;
  if (cws.size() < required || ys.size() < kOutputsPerKey ||
      cws.data() == nullptr || ys.data() == nullptr ||
      kOutputsPerKey > std::numeric_limits<size_t>::max() / sizeof(int4))
    return cudaErrorInvalidValue;
  gpu::detail::PackedHalfTreeSplitEvalAllKernel<Scheme, LaunchShape::kFrontierDepth,
      LaunchShape::kBlockRootDepth>
      <<<LaunchShape::kBlocksPerKey, LaunchShape::kThreadsPerBlock, 0, stream>>>(
          b, s0, cws.data(), stride, ys.data(), instance_seed_, hash_);
  return cudaGetLastError();
#else
  static_assert(sizeof(Hash) == 0, "gpu launch requires a CUDA compiler");
  return cudaErrorNotSupported;
#endif
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
template <int frontier_depth, int threads_per_block>
__host__ inline cudaError_t PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalAllGpuBatch(
    bool b, cuda::std::span<const int4> s0s, cuda::std::span<const int4> cws,
    cuda::std::span<int4> ys, cudaStream_t stream, size_t stride) const
  requires(kDepth < sizeof(size_t) * 8) {
#if defined(__CUDACC__)
  using LaunchShape = GpuLaunchShape<frontier_depth, threads_per_block>;
  using Scheme = PackedHalfTreeDpf<in_bits, out_bits, Hash, In>;
  constexpr size_t kOutputsPerKey = size_t{1} << kDepth;
  constexpr size_t kMaxElements = std::numeric_limits<size_t>::max() / sizeof(int4);
  constexpr size_t kMaxAddressableSize = std::numeric_limits<size_t>::max();
  const size_t key_count = s0s.size();
  if (key_count == 0) return cudaSuccess;
  // The output guard bounds key_count * kOutputsPerKey * sizeof(int4) so the
  // product in the size check below cannot wrap.
  if (key_count > kMaxElements ||
      key_count > kMaxAddressableSize / kOutputsPerKey / sizeof(int4) ||
      (stride != 0 && (stride < key_count || stride > kMaxElements)))
    return cudaErrorInvalidValue;
  // A stride of 0 resolves to the key count, which is the dense level-major layout.
  const size_t resolved_stride = stride == 0 ? key_count : stride;
  if (resolved_stride > (kMaxElements - key_count) / static_cast<size_t>(kDepth))
    return cudaErrorInvalidValue;
  const size_t required = static_cast<size_t>(kDepth) * resolved_stride + key_count;
  if (cws.size() < required || ys.size() < key_count * kOutputsPerKey ||
      s0s.data() == nullptr || cws.data() == nullptr || ys.data() == nullptr)
    return cudaErrorInvalidValue;
  if (key_count > static_cast<size_t>(std::numeric_limits<int>::max()) / LaunchShape::kBlocksPerKey)
    return cudaErrorInvalidConfiguration;
  gpu::detail::PackedHalfTreeSplitEvalAllBatchKernel<Scheme, LaunchShape::kFrontierDepth,
      LaunchShape::kBlockRootDepth>
      <<<key_count * LaunchShape::kBlocksPerKey, LaunchShape::kThreadsPerBlock, 0, stream>>>(
          b, s0s.data(), cws.data(), resolved_stride, ys.data(), instance_seed_, hash_);
  return cudaGetLastError();
#else
  static_assert(sizeof(Hash) == 0, "gpu launch requires a CUDA compiler");
  return cudaErrorNotSupported;
#endif
}
 
#if FSS_ENABLE_V1_COMPAT
 
// Deprecated v1.3.0 aggregate-Key entries. Every wrapper copies between the
// caller's Key layout and a compile-time-sized split-component scratch buffer,
// then runs the current public method at stride 1. The old layout kept the
// packed-output correction word in the caller's separate fcw.
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::Gen(
    compat::detail::FcwKeyRecord<PackedHalfTreeDpf, Cw, kDepth> &k0,
    compat::detail::FcwKeyRecord<PackedHalfTreeDpf, Cw, kDepth> &k1,
    cuda::std::span<const int4, 2> s0s, In a, uint64_t beta) {
  int4 seed[kDepth + 1];
  Gen(cuda::std::span<int4>(seed, kDepth + 1), s0s, a, beta);
  for (int i = 0; i < kDepth; ++i) {
    const Cw layer = {seed[i]};
    k0.cws[i] = layer;
    k1.cws[i] = layer;
  }
  k0.fcw = seed[kDepth];
  k1.fcw = seed[kDepth];
  k0.s0 = s0s[0];
  k1.s0 = s0s[1];
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
__host__ __device__ inline int4 PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::Eval(
    bool b, const compat::detail::FcwKeyRecord<PackedHalfTreeDpf, Cw, kDepth> &key, In x) {
  int4 seed[kDepth + 1];
  for (int i = 0; i < kDepth; ++i) seed[i] = key.cws[i].s;
  seed[kDepth] = key.fcw;
  return Eval(b, key.s0, cuda::std::span<const int4>(seed, kDepth + 1), x);
}
 
template <int in_bits, int out_bits, typename Hash, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      CcrHashable<Hash> && out_bits >= 1 && out_bits <= 64 && (out_bits & (out_bits - 1)) == 0)
inline void PackedHalfTreeDpf<in_bits, out_bits, Hash, In>::EvalAll(
    bool b, const compat::detail::FcwKeyRecord<PackedHalfTreeDpf, Cw, kDepth> &key,
    cuda::std::span<int4> ys) {
  int4 seed[kDepth + 1];
  for (int i = 0; i < kDepth; ++i) seed[i] = key.cws[i].s;
  seed[kDepth] = key.fcw;
  EvalAll(b, key.s0, cuda::std::span<const int4>(seed, kDepth + 1), ys.data());
}
 
#endif  // FSS_ENABLE_V1_COMPAT
 
}  // namespace fss