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>
#ifdef _OPENMP
#include <omp.h>
#endif
#include <vector>
#include <cuda/std/array>
#include <fss/group.cuh>
#include <fss/prg.cuh>
#include <fss/util.cuh>
#include <cuda/std/span>
#include <limits>
#include <fss/constraints.cuh>
#include <fss/compat.cuh>
 
#include <cstdint>
#include <stdexcept>
 
namespace fss {
 
namespace detail {
class SchemeAccess;
}  // namespace detail
 
namespace detail {
 
// Split correction-word view shared by the host and device split paths. Layer i is
// read from cws[i*stride] and the right-control correction is bit i of
// packed_tr. The terminal is rebuilt as a plain Cw with no control bit, so bits of
// packed_tr at or above in_bits are never read; the shift by in_bits would be
// undefined when in_bits equals the width of In. The layer LSB keeps the
// left-control correction exactly as the Cw layout does. operator[] recovers the Cw
// the shared algorithms consume, so their bodies stay identical on the Cw and split
// paths.
template <typename Cw, typename In, int in_bits>
struct DpfCwsView {
  const int4 *s;
  In packed_tr;
  size_t stride;
 
  __host__ __device__ inline Cw operator[](int i) const {
    const int4 s_cw = s[static_cast<size_t>(i) * stride];
    if (i >= in_bits) return Cw{s_cw, false};
    return Cw{s_cw, ((packed_tr >> i) & 1) != 0};
  }
};
 
}  // namespace detail
 
template <int in_bits, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> && Prgable<Prg, 2>)
class GrottoDcf;

template <int in_bits, typename Group, typename Prg, typename In = uint>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
class Dpf {
public:
  struct __align__(32) Cw {
    int4 s;  
    bool tr;  
  };
  // For only 1 and aligned memory access on GPU
  static_assert(sizeof(Cw) == 32);

  struct Key {
    Cw cws[in_bits + 1];  
    int4 s0;  
#if FSS_ENABLE_V1_COMPAT
    [[deprecated("the par_depth template parameter is gone; pass par_depth to EvalAll() instead")]]
    static constexpr int kParDepth = -1;
#endif  // FSS_ENABLE_V1_COMPAT
  };

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

  __host__ __device__ explicit Dpf(Prg prg);

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

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

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

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

  template <int frontier_depth = 16, int threads_per_block = 256>
  __host__ cudaError_t EvalRangeGpu(bool b, In first, In last, int4 s0,
      cuda::std::span<const int4> cws, In tr, cuda::std::span<int4> ys,
      cudaStream_t stream = 0, size_t stride = 1) const;

  template <int frontier_depth = 16, int threads_per_block = 256>
  __host__ cudaError_t EvalAllGpu(bool b, int4 s0, cuda::std::span<const int4> cws,
      In tr, cuda::std::span<int4> ys, cudaStream_t stream = 0, size_t stride = 1) const
    requires(in_bits < 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<const In> tr,
      cuda::std::span<int4> ys, cudaStream_t stream = 0, size_t stride = 0) const
    requires(in_bits < sizeof(size_t) * 8);
 
#if FSS_ENABLE_V1_COMPAT
  // Deprecated v1.3.0 aggregate-Cw compatibility entries; see fss/compat.cuh.
  // Each converts through a compile-time-sized scratch buffer and the current
  // split-component method above at stride 1. cws has in_bits + 1 entries.

  [[deprecated("use the split-component Gen(cuda::std::span<int4>, In&, cuda::std::span<const int4, 2>, In, int4, size_t) overload")]]
  __host__ __device__ void Gen(Cw cws[], const int4 s0s[2], In a, int4 b_buf);

  [[deprecated("use the split-component Eval(bool, int4, cuda::std::span<const int4>, In, In, size_t) overload")]]
  __host__ __device__ int4 Eval(bool b, int4 s0, const Cw cws[], In x);

  [[deprecated("use the split-component EvalAll(bool, int4, cuda::std::span<const int4>, In, int4[], int, size_t) overload")]]
  void EvalAll(bool b, int4 s0, const Cw cws[], int4 ys[]);
 
  // Deprecated v1.3.0 aggregate-Key entries. They keep the current public Key
  // type, which the modern bucket APIs still use, and convert through the same
  // compile-time scratch buffer as the Cw entries above. EvalAll() runs
  // sequentially, as the current split-component default does; the old
  // automatic thread-count split is not restored.

  [[deprecated("use the split-component Gen(cuda::std::span<int4>, In&, cuda::std::span<const int4, 2>, In, int4, size_t) overload")]]
  __host__ __device__ void Gen(Key &k0, Key &k1, cuda::std::span<const int4, 2> s0s,
      In a, int4 b_buf);

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

  [[deprecated("use the split-component EvalAll(bool, int4, cuda::std::span<const int4>, In, int4[], int, size_t) overload")]]
  void EvalAll(bool b, const Key &key, cuda::std::span<int4> ys);
#endif  // FSS_ENABLE_V1_COMPAT
 
private:
  friend class detail::SchemeAccess;
 
  Prg prg_;
 
  template <int frontier_depth, int threads_per_block>
  struct GpuLaunchShape;
 
  // EvalAll() shape: a scalar prefix of kBatchPrefix levels grows the frontier to
  // 8 nodes, detail::GenBatch() expands 8 parents per level (pipelined through
  // Gen8() when the PRG offers a batched entry), and a kBatchTail-level tail
  // converts each 8-subtree group of 8 << kBatchTail leaves in local memory.
  static constexpr int kBatchPrefix = 3;
  static constexpr int kBatchTail = 3;
 
  friend class GrottoDcf<in_bits, Prg, In>;
 
  // Packed states store the control bit in their least significant bit.
  // correction retains the left control bit in its LSB; right_control is Cw::tr.
  // A backend reference lets CPU and GPU callers reuse the same calculation
  // without reconstructing a Scheme or copying its backend at every node.
  __host__ __device__ static void ExpandNode(Prg &prg, int4 state, int4 correction,
      bool right_control, int4 &left, int4 &right);
 
  __host__ __device__ static int4 PathStep(Prg &prg, int4 state, int4 correction,
      bool right_control, bool take_right);
 
  // Correction is a raw int4, or the Group already converted by the CPU batch
  // tail. Raw correction values are converted only when the leaf control bit is set.
  template <typename Correction>
  __host__ __device__ static int4 ConvertLeaf(bool party, int4 state, Correction &correction);
 
  // Split correction-word view over the caller's arrays. The host methods and the
  // split GPU kernels share this one view; see detail::DpfCwsView.
  using CwsView = detail::DpfCwsView<Cw, In, in_bits>;
 
  // Writes Cw slots in the legacy AoS layout, covering each full 32B slot.
  struct AosWriter {
    Cw *cws;
 
    __host__ __device__ void Write(int i, int4 s_cw, bool tr_cw) const {
      // Write both 16B halves so the full 32B slot is covered by unmasked stores.
      // {s_cw, tr_cw} alone leaves the padding undefined: NVCC emits a masked
      // 256-bit store, and the partial-sector write forces a write-allocate
      // round trip per level on GPU (measured ~2x Gen time).
      int4 *cw_slot = reinterpret_cast<int4 *>(&cws[i]);
      cw_slot[0] = s_cw;
      cw_slot[1] = make_int4(tr_cw, 0, 0, 0);
    }
 
    __host__ __device__ void Finish(int4 terminal) const {
      cws[in_bits] = {terminal, false};
    }
  };
 
  // Writes strided seed components and accumulates packed right-control bits locally.
  struct CwsWriter {
    int4 *cws;
    size_t stride;
    In &tr;
    In packed_tr = 0;
 
    __host__ __device__ void Write(int i, int4 s_cw, bool tr_cw) {
      cws[static_cast<size_t>(i) * stride] = s_cw;
      packed_tr |= static_cast<In>(tr_cw) << i;
    }
 
    __host__ __device__ void Finish(int4 terminal) {
      cws[static_cast<size_t>(in_bits) * stride] = terminal;
      tr = packed_tr;
    }
  };
 
  // Shared correction 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, int4 b_buf);
 
  // Walks x's path from the root and returns the leaf seed with its control bit in the LSB.
  // 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);
 
  // Shared scalar evaluation over the Cw array or the split view.
  template <typename Cws>
  __host__ __device__ int4 EvalCws(bool b, int4 s0, Cws cws, In x);
 
  // Shared interval evaluation over the Cw array or the split view.
  template <typename Cws>
  void EvalRangeCws(bool b, In first, In last, int4 s0, Cws cws,
      cuda::std::span<int4> ys, int par_depth);
 
  // Shared full-domain evaluation over the Cw array or the split view.
  template <typename Cws>
  void EvalAllCws(bool b, int4 s0, Cws cws, int4 ys[], int par_depth);
 
  // Full-domain batch GPU launch over key-major rows of Cw, the storage
  // contract of DMPF bucket keys. Reached through detail::SchemeAccess only;
  // public callers use the level-major split batch API.
  template <int frontier_depth, int threads_per_block>
  __host__ cudaError_t EvalAllGpuBatchCw(bool b, cuda::std::span<const int4> s0s,
      cuda::std::span<const Cw> cws, cuda::std::span<int4> ys, cudaStream_t stream) const;
 
  // Closed node bounds avoid representing an unaddressable full-domain end.
  template <typename Cws>
  void EvalRangeSerial(bool b, In first, In last, int4 s0, Cws cws, int4 ys[]);
 
  // GCC 12.4.0 can discard these writes before the OpenMP walk after IPA
  // effect analysis. Keep this boundary opaque on GCC 12 releases before 12.5.
  // GCC PR115033 describes a related closure escape-analysis defect.
  template <typename Cws>
#if defined(__GNUC__) && !defined(__clang__) && __GNUC__ == 12 && __GNUC_MINOR__ <= 4
  __attribute__((noipa))
#endif
  std::vector<int4> BuildParallelFrontier(int4 st, Cws cws, int split_depth);
 
  // Expands w scalar frontier nodes (t in LSB) by one level with cws[lv].
  // The children overwrite the frontier in place: iteration is in reverse so
  // the writes at 2j, 2j+1 never clobber nodes still to be read.
  template <typename Cws>
  void ExpandFrontier(int4 frontier[], size_t w, Cws cws, int lv);
 
  // Applies the level correction to one expanded child branchlessly.
  // m is 0 or -1 from the parent control bit t, and t_cw is the cw control bit
  // masked by t. The cw seed is clamped, so the XOR leaves the child LSB free
  // for the packed control bit.
  __host__ __device__ static int4 CorrectChild(int4 child, int4 s_cw, int m, int t_cw);
 
  // Full-domain evaluation of the subtree rooted at st (t in LSB) at level i through an
  // 8-wide batched frontier: a kBatchPrefix-level scalar prefix, batched expansions, and
  // a kBatchTail-level tail converted in local memory.
  // Requires in_bits - i >= kBatchPrefix + kBatchTail.
  template <typename Cws>
  void EvalBatched8(bool b, int4 st, Cws cws, int4 ys[], size_t l, int i);
 
  // Evaluates a group of 8 consecutive frontier nodes (t in LSB) at level i, covering the
  // leaves [l, l + 8 << (in_bits - i)) in order. Requires in_bits - i >= kBatchTail.
  template <typename Cws>
  void EvalBatched8Group(bool b, const int4 (&frontier)[8], Cws cws, int4 ys[], size_t l, int i);
};
 
// -----------------------------------------------------------------------------
 
namespace gpu::detail {
 
template <typename Scheme>
struct DpfGpuSchemeTraits;
 
#if defined(__CUDACC__)
 
// The single-key and batch kernels share the same block-local frontier and tail.
// Each node keeps the former full-domain helper's by-value backend lifetime.
// Cws is a Cw pointer on the array path and a kernel-boundary split view
// on the split path, so the traversal body does not depend on the layout.
template <typename Scheme, int frontier_depth, int block_root_depth, bool range = false, typename Cws>
__device__ __forceinline__ void DpfEvalAllBlock(bool b, int4 s0, Cws cws,
    int4 *ys, typename DpfGpuSchemeTraits<Scheme>::Backend &prg, int subtree_block_index, int64_t thread_index, int4 *shared_frontier,
    __uint128_t first = 0, int subtree_depth = DpfGpuSchemeTraits<Scheme>::kInputBits,
    int range_frontier = frontier_depth, int range_block_root = block_root_depth) {
  constexpr int kInputBits = DpfGpuSchemeTraits<Scheme>::kInputBits;
  constexpr int kTailCapacity = range ? 8 : kInputBits - frontier_depth;
  const int start_depth = range ? kInputBits - subtree_depth : 0;
  const int frontier = range ? start_depth + range_frontier : frontier_depth;
  const int root_depth = range ? start_depth + range_block_root : block_root_depth;
  const int tail_depth = kInputBits - frontier;
  const int block_subtree_depth = frontier - root_depth;
  // Phase 0: block root walk (top block_root_depth levels) by thread 0.
  if (threadIdx.x == 0) {
    int4 st = util::SetLsb(s0, b);
#pragma unroll
    for (int i = 0; i < root_depth; ++i) {
      bool x_bit = i < start_depth
          ? ((first >> (kInputBits - 1 - i)) & 1) != 0
          : ((subtree_block_index >> (root_depth - 1 - i)) & 1) != 0;
      {
        // Preserve the old by-value path helper's copy and destruction boundary.
        auto node_prg = prg;
        st = fss::detail::SchemeAccess::DpfPathStep<Scheme>(node_prg, st, cws[i].s, cws[i].tr, x_bit);
      }
    }
    shared_frontier[0] = st;
  }
  __syncthreads();
 
  // Phase 1: breadth-parallel expansion of the block's subtree, levels block_root_depth..frontier_depth.
#pragma unroll
  for (int i = 0; i < block_subtree_depth; ++i) {
    const bool active = threadIdx.x < (1 << (i + 1));
    int4 n = {0, 0, 0, 0};
    if (active) n = shared_frontier[threadIdx.x >> 1];
    // All parent reads must finish before any child overwrites the frontier.
    __syncthreads();
    if (active) {
      int4 left, right;
      {
        // Preserve the old by-value expansion helper's copy and destruction boundary.
        auto node_prg = prg;
        fss::detail::SchemeAccess::DpfExpandNode<Scheme>(node_prg, n, cws[root_depth + i].s, cws[root_depth + i].tr, left, right);
      }
      shared_frontier[threadIdx.x] = (threadIdx.x & 1) ? right : left;
    }
    __syncthreads();
  }
 
  if constexpr (range) {
    int4 pending[8];
    int4 node = shared_frontier[threadIdx.x];
    int level = frontier;
    const int count = 1 << tail_depth;
    const size_t offset = static_cast<size_t>(thread_index) << tail_depth;
    for (int index = 0; index < count; ++index) {
      while (level < kInputBits) {
        int4 left, right;
        {
          auto node_prg = prg;
          fss::detail::SchemeAccess::DpfExpandNode<Scheme>(node_prg, node,
              cws[level].s, cws[level].tr, left, right);
        }
        pending[level - frontier] = right;
        node = left;
        ++level;
      }
      ys[offset + index] = fss::detail::SchemeAccess::DpfConvertLeaf<Scheme>(
          b, node, cws[kInputBits].s);
      if (index + 1 < count) {
        int next = index + 1;
        level = kInputBits - 1;
        while ((next & 1) == 0) {
          --level;
          next >>= 1;
        }
        node = pending[level - frontier];
        ++level;
      }
    }
  } else {
    // Phase 2: expand the subtree below depth frontier_depth, then convert to leaves.
    int4 nodes[1 << kTailCapacity];
    nodes[0] = shared_frontier[threadIdx.x];
    int count = 1;
#pragma unroll
    for (int lev = frontier; lev < kInputBits; ++lev) {
#pragma unroll
      for (int k = count - 1; k >= 0; --k) {
        int4 left, right;
        {
          // Preserve the old by-value expansion helper's copy and destruction boundary.
          auto node_prg = prg;
          fss::detail::SchemeAccess::DpfExpandNode<Scheme>(node_prg, nodes[k], cws[lev].s, cws[lev].tr, left, right);
        }
        nodes[2 * k] = left;
        nodes[2 * k + 1] = right;
      }
      count *= 2;
    }
 
    int4 v_cw_np1 = cws[kInputBits].s;
    const size_t output_offset = static_cast<size_t>(thread_index) << tail_depth;
#pragma unroll
    for (int k = 0; k < count; ++k) {
      ys[output_offset + k] = fss::detail::SchemeAccess::DpfConvertLeaf<Scheme>(b, nodes[k], v_cw_np1);
    }
  }
}
 
template <typename Scheme, int frontier_depth, int block_root_depth>
__global__ void DpfEvalAllBatchKernel(bool b, const int4 *s0s, const void *cws_buf,
    int4 *ys, typename DpfGpuSchemeTraits<Scheme>::Backend prg) {
  // The caller passes the raw int4 correction records; the array layout is one
  // Cw per record, so the typed view is a plain cast.
  const auto *cws = static_cast<const typename Scheme::Cw *>(cws_buf);
  constexpr int kInputBits = DpfGpuSchemeTraits<Scheme>::kInputBits;
  constexpr int kBlocksPerKey = 1 << block_root_depth;
  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 auto *key_cws = cws + static_cast<size_t>(key_index) * (kInputBits + 1);
  int4 *key_outputs = ys + static_cast<size_t>(key_index) * (size_t{1} << kInputBits);
  __shared__ int4 shared_frontier[1 << (frontier_depth - block_root_depth)];
  DpfEvalAllBlock<Scheme, frontier_depth, block_root_depth>(b, s0s[key_index], key_cws, key_outputs, prg, subtree_block_index, thread_index, shared_frontier);
}
 
// Split-layout counterparts of the three kernels above. The correction buffers
// arrive as raw records and the boundary view reconstructs the per-layer Cw
// entries, so the shared traversal stays identical on both paths.
template <typename Scheme, int threads_per_block>
__global__ void DpfSplitEvalRangeKernel(bool b, int4 s0,
    const int4 *cws, typename DpfGpuSchemeTraits<Scheme>::Input tr,
    size_t stride, int4 *ys, __uint128_t first,
    int depth, int frontier, int root, typename DpfGpuSchemeTraits<Scheme>::Backend prg) {
  using Cws = fss::detail::DpfCwsView<typename Scheme::Cw,
      typename DpfGpuSchemeTraits<Scheme>::Input,
      DpfGpuSchemeTraits<Scheme>::kInputBits>;
  __shared__ int4 states[threads_per_block];
  const int64_t thread_index = static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
  DpfEvalAllBlock<Scheme, 0, 0, true>(b, s0, Cws{cws, tr, stride}, ys, prg,
      blockIdx.x, thread_index, states, first, depth, frontier, root);
}
 
template <typename Scheme, int frontier_depth, int block_root_depth>
__global__ void DpfSplitEvalAllKernel(bool b, int4 s0, const int4 *cws,
    typename DpfGpuSchemeTraits<Scheme>::Input tr, size_t stride,
    int4 *ys, typename DpfGpuSchemeTraits<Scheme>::Backend prg) {
  using Cws = fss::detail::DpfCwsView<typename Scheme::Cw,
      typename DpfGpuSchemeTraits<Scheme>::Input,
      DpfGpuSchemeTraits<Scheme>::kInputBits>;
  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)];
  DpfEvalAllBlock<Scheme, frontier_depth, block_root_depth>(b, s0, Cws{cws, tr, stride},
      ys, prg, blockIdx.x, thread_index, shared_frontier);
}
 
template <typename Scheme, int frontier_depth, int block_root_depth>
__global__ void DpfSplitEvalAllBatchKernel(bool b, const int4 *s0s, const int4 *cws,
    const typename DpfGpuSchemeTraits<Scheme>::Input *trs, size_t stride,
    int4 *ys, typename DpfGpuSchemeTraits<Scheme>::Backend prg) {
  using Cws = fss::detail::DpfCwsView<typename Scheme::Cw,
      typename DpfGpuSchemeTraits<Scheme>::Input,
      DpfGpuSchemeTraits<Scheme>::kInputBits>;
  constexpr int kInputBits = DpfGpuSchemeTraits<Scheme>::kInputBits;
  constexpr int kBlocksPerKey = 1 << block_root_depth;
  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) * (size_t{1} << kInputBits);
  __shared__ int4 shared_frontier[1 << (frontier_depth - block_root_depth)];
  DpfEvalAllBlock<Scheme, frontier_depth, block_root_depth>(b, s0s[key_index],
      Cws{cws + key_index, trs[key_index], stride}, key_outputs, prg,
      subtree_block_index, thread_index, shared_frontier);
}
 
// Device-key counterparts of the split kernels above: the party seed and the
// packed control live in caller device storage and are loaded here, so no host
// value of the key material exists. They feed the same shared traversal.
template <typename Scheme, int frontier_depth, int block_root_depth>
__global__ void DpfSplitEvalAllDeviceKeyKernel(bool b, const int4 *s0,
    const int4 *cws, const typename DpfGpuSchemeTraits<Scheme>::Input *tr,
    size_t stride, int4 *ys, typename DpfGpuSchemeTraits<Scheme>::Backend prg) {
  using Cws = fss::detail::DpfCwsView<typename Scheme::Cw,
      typename DpfGpuSchemeTraits<Scheme>::Input,
      DpfGpuSchemeTraits<Scheme>::kInputBits>;
  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)];
  DpfEvalAllBlock<Scheme, frontier_depth, block_root_depth>(b, s0[0],
      Cws{cws, tr[0], stride}, ys, prg, blockIdx.x, thread_index, shared_frontier);
}
 
template <typename Scheme, int threads_per_block>
__global__ void DpfSplitEvalRangeDeviceKeyKernel(bool b, const int4 *s0,
    const int4 *cws, const typename DpfGpuSchemeTraits<Scheme>::Input *tr,
    size_t stride, int4 *ys, __uint128_t first,
    int depth, int frontier, int root, typename DpfGpuSchemeTraits<Scheme>::Backend prg) {
  using Cws = fss::detail::DpfCwsView<typename Scheme::Cw,
      typename DpfGpuSchemeTraits<Scheme>::Input,
      DpfGpuSchemeTraits<Scheme>::kInputBits>;
  __shared__ int4 states[threads_per_block];
  const int64_t thread_index = static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x;
  DpfEvalAllBlock<Scheme, 0, 0, true>(b, s0[0], Cws{cws, tr[0], stride}, ys, prg,
      blockIdx.x, thread_index, states, first, depth, frontier, root);
}
 
// Key-major batch counterpart: key k's components occupy records
// [k*key_records, k*key_records + in_bits] and its control record follows, so
// the per-key view uses a unit stride.
template <typename Scheme, int frontier_depth, int block_root_depth>
__global__ void DpfSplitEvalAllBatchDeviceKeyKernel(bool b, const int4 *s0s,
    const int4 *cws, size_t key_records, size_t control_record,
    size_t key_count, int4 *ys, typename DpfGpuSchemeTraits<Scheme>::Backend prg) {
  using In = typename DpfGpuSchemeTraits<Scheme>::Input;
  using Cws = fss::detail::DpfCwsView<typename Scheme::Cw, In,
      DpfGpuSchemeTraits<Scheme>::kInputBits>;
  constexpr int kInputBits = DpfGpuSchemeTraits<Scheme>::kInputBits;
  constexpr int kBlocksPerKey = 1 << block_root_depth;
  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;
  const In *tr = reinterpret_cast<const In *>(key + control_record);
  int4 *key_outputs = ys + static_cast<size_t>(key_index) * (size_t{1} << kInputBits);
  __shared__ int4 shared_frontier[1 << (frontier_depth - block_root_depth)];
  DpfEvalAllBlock<Scheme, frontier_depth, block_root_depth>(b, s0s[key_index],
      Cws{key, tr[0], 1}, key_outputs, prg, subtree_block_index, thread_index, shared_frontier);
}
 
// Host launchers for the device-key kernels. The Prg is passed by value from a
// host mirror of the device scheme, matching the public launch methods.
template <int in_bits, typename Group, typename Prg, typename In, int frontier_depth = 16,
    int threads_per_block = 256>
__host__ inline cudaError_t DpfEvalAllGpuDeviceKey(bool b, const int4 *s0,
    const int4 *cws, size_t cws_count, const In *tr, size_t stride, int4 *ys, size_t ys_count,
    const Prg &prg, cudaStream_t stream) {
  constexpr size_t kMaxElements = std::numeric_limits<size_t>::max() / sizeof(int4);
  constexpr size_t kOutputsPerKey = size_t{1} << in_bits;
  constexpr int kFrontierDepth = frontier_depth < in_bits ? frontier_depth : in_bits;
  if (stride == 0 || stride > (kMaxElements - 1) / static_cast<size_t>(in_bits) ||
      cws_count < static_cast<size_t>(in_bits) * stride + 1 || ys_count < kOutputsPerKey ||
      s0 == nullptr || cws == nullptr || tr == nullptr || ys == nullptr ||
      kOutputsPerKey > kMaxElements)
    return cudaErrorInvalidValue;
  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(in_bits - 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;
  using Scheme = Dpf<in_bits, Group, Prg, In>;
  DpfSplitEvalAllDeviceKeyKernel<Scheme, kFrontierDepth, kBlockRootDepth>
      <<<size_t{1} << kBlockRootDepth, 1 << kThreadsPerBlockLog2, 0, stream>>>(
          b, s0, cws, tr, stride, ys, prg);
  return cudaGetLastError();
}
 
template <int in_bits, typename Group, typename Prg, typename In, int frontier_depth = 16,
    int threads_per_block = 256>
__host__ inline cudaError_t DpfEvalRangeGpuDeviceKey(bool b, In first, In last,
    const int4 *s0, const int4 *cws, size_t cws_count, const In *tr, size_t stride,
    int4 *ys, size_t ys_count, const Prg &prg, cudaStream_t stream) {
  constexpr size_t kMaxElements = std::numeric_limits<size_t>::max() / sizeof(int4);
  if (stride == 0 || stride > (kMaxElements - 1) / static_cast<size_t>(in_bits) ||
      cws_count < static_cast<size_t>(in_bits) * stride + 1 || cws == nullptr ||
      s0 == nullptr || tr == nullptr || ys == nullptr) return cudaErrorInvalidValue;
  using Scheme = Dpf<in_bits, Group, Prg, In>;
  return fss::detail::SplitRangeChunks<in_bits, frontier_depth, threads_per_block>(
      first, last, ys_count,
      [&](In chunk_first, int depth, size_t offset, int frontier, int thread_bits, int root,
          size_t blocks) {
        DpfSplitEvalRangeDeviceKeyKernel<Scheme, threads_per_block>
            <<<blocks, 1 << thread_bits, 0, stream>>>(b, s0, cws, tr, stride,
                ys + offset, chunk_first, depth, frontier, root, prg);
        return cudaGetLastError();
      });
}
 
template <int in_bits, typename Group, typename Prg, typename In, int frontier_depth = 16,
    int threads_per_block = 256>
__host__ inline cudaError_t DpfEvalAllGpuBatchDeviceKey(bool b, const int4 *s0s,
    size_t key_count, const int4 *cws, size_t cws_count, size_t key_records,
    size_t control_record, int4 *ys, size_t ys_count, const Prg &prg, cudaStream_t stream) {
  constexpr size_t kMaxElements = std::numeric_limits<size_t>::max() / sizeof(int4);
  constexpr size_t kOutputsPerKey = size_t{1} << in_bits;
  constexpr size_t kMaxAddressableSize = std::numeric_limits<size_t>::max();
  constexpr int kFrontierDepth = frontier_depth < in_bits ? frontier_depth : in_bits;
  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(in_bits - 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 > 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 ||
      kOutputsPerKey > kMaxElements) return cudaErrorInvalidValue;
  if (key_count > static_cast<size_t>(std::numeric_limits<int>::max()) / kBlocksPerKey)
    return cudaErrorInvalidConfiguration;
  using Scheme = Dpf<in_bits, Group, Prg, In>;
  DpfSplitEvalAllBatchDeviceKeyKernel<Scheme, kFrontierDepth, kBlockRootDepth>
      <<<key_count * kBlocksPerKey, 1 << kThreadsPerBlockLog2, 0, stream>>>(
          b, s0s, cws, key_records, control_record, key_count, ys, prg);
  return cudaGetLastError();
}
 
// Correction words arrive in the dense level-major split layout: layer i of
// key k is cws[i * key_count + k], the terminal layer is
// cws[in_bits * key_count + k], and tr[k] holds key k's packed right-control
// bits, bit i for layer i.
template <typename Scheme, int unroll, bool b, typename Control, typename Query>
__global__ void DpfCoreEvalPointKernel(int4 *ys, const int4 *seeds, const int4 *cws,
    const Control *tr, const int4 *out_cw,
    const Query *xs, int key_count,
    typename DpfGpuSchemeTraits<Scheme>::Backend prg) {
  using In = typename DpfGpuSchemeTraits<Scheme>::Input;
  constexpr int kInputBits = DpfGpuSchemeTraits<Scheme>::kInputBits;
  static_assert(kInputBits <= sizeof(In) * 8, "In too narrow for kInputBits");
  int thread_index = blockIdx.x * blockDim.x + threadIdx.x;
  if (thread_index >= key_count) return;
  int4 state = util::SetLsb(seeds[thread_index], b);
  In x = static_cast<In>(xs[thread_index]);
  Control ex = tr[thread_index];
#pragma unroll unroll
  for (int i = 0; i < kInputBits; ++i) {
    bool take_right = ((x >> (kInputBits - 1 - i)) & 1) != 0;
    bool right_control = ((ex >> i) & 1) != 0;
    state = fss::detail::SchemeAccess::DpfPathStep<Scheme>(
        prg, state, cws[static_cast<size_t>(i) * key_count + thread_index], right_control, take_right);
  }
  ys[thread_index] = fss::detail::SchemeAccess::DpfConvertLeaf<Scheme>(b, state, out_cw[thread_index]);
}
 
template <typename Scheme, typename Control, typename Query>
struct DpfPointGpuLauncher;
 
template <int n, typename Group, typename Prg, typename In, typename Control, typename Query>
struct DpfPointGpuLauncher<Dpf<n, Group, Prg, In>, Control, Query> {
  template <int unroll>
  static void Launch(bool b, const int4 *seeds, const int4 *cws, const Control *tr,
      const int4 *out_cw, const Query *xs, int4 *ys, int key_count, Prg prg,
      cudaStream_t stream) {
    using Scheme = Dpf<n, Group, Prg, In>;
    constexpr int kThreadsPerBlock = 256;
    int block_count = (key_count + kThreadsPerBlock - 1) / kThreadsPerBlock;
    if (b) {
      DpfCoreEvalPointKernel<Scheme, unroll, true, Control, Query>
          <<<block_count, kThreadsPerBlock, 0, stream>>>(ys, seeds, cws, tr, out_cw, xs, key_count, prg);
    } else {
      DpfCoreEvalPointKernel<Scheme, unroll, false, Control, Query>
          <<<block_count, kThreadsPerBlock, 0, stream>>>(ys, seeds, cws, tr, out_cw, xs, key_count, prg);
    }
  }
};
 
template <int unroll, typename Scheme, typename Control, typename Query>
void DpfCoreEvalPointLaunch(bool b, const int4 *seeds, const int4 *cws, const Control *tr,
    const int4 *out_cw, const Query *xs, int4 *ys, int key_count,
    typename DpfGpuSchemeTraits<Scheme>::Backend prg, cudaStream_t stream) {
  if (key_count == 0) return;
  DpfPointGpuLauncher<Scheme, Control, Query>::template Launch<unroll>(
      b, seeds, cws, tr, out_cw, xs, ys, key_count, prg, stream);
}

template <int unroll = 4, int in_bits, typename Group, typename Prg, typename In, typename Query = In>
  requires QueryFor<Query, in_bits>
void DpfEvalPointGpu(bool b, cuda::std::span<const int4> seeds, cuda::std::span<const int4> cws,
    cuda::std::span<const In> tr, cuda::std::span<const Query> xs, int4 *ys, int nkeys,
    const Dpf<in_bits, Group, Prg, In> &dpf, cudaStream_t stream = 0) {
  if (nkeys < 0)
    throw std::invalid_argument("DpfEvalPointGpu: 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>(in_bits) + 1;
  if (count > SIZE_MAX / kLevels || cws.size() < count * kLevels)
    throw std::invalid_argument("DpfEvalPointGpu: cws must hold (in_bits+1)*nkeys elements");
  if (seeds.size() < count || tr.size() < count || xs.size() < count || ys == nullptr)
    throw std::invalid_argument("DpfEvalPointGpu: seeds, tr, xs, and ys must cover nkeys entries");
  const int4 *out_cw = cws.data() + static_cast<size_t>(in_bits) * static_cast<size_t>(nkeys);
  DpfCoreEvalPointLaunch<unroll, Dpf<in_bits, Group, Prg, In>, In, Query>(b, seeds.data(), cws.data(), tr.data(),
      out_cw, xs.data(), ys, nkeys, fss::detail::SchemeAccess::GetPrg(dpf), stream);
}
 
#endif  // defined(__CUDACC__)
 
}  // namespace gpu::detail
 
namespace gpu::detail {
 
template <int n, typename Group, typename Prg, typename In>
struct DpfGpuSchemeTraits<Dpf<n, Group, Prg, In>> {
  using Backend = Prg;
  using Input = In;
  static constexpr int kInputBits = n;
};
 
}  // namespace gpu::detail
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline Dpf<in_bits, Group, Prg, In>::Dpf(Prg prg)
    : prg_(prg) {}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline void Dpf<in_bits, Group, Prg, In>::Gen(
    cuda::std::span<int4> cws, In &tr, cuda::std::span<const int4, 2> s0s,
    In a, int4 b_buf, size_t stride) {
  util::RequireArg(stride > 0, "Dpf::Gen: stride must be positive");
  constexpr size_t kMaxElements = SIZE_MAX / sizeof(int4);
  util::RequireArg(stride <= (kMaxElements - 1) / static_cast<size_t>(in_bits),
      "Dpf::Gen: stride exceeds the correction array limit");
  const size_t required = static_cast<size_t>(in_bits) * stride + 1;
  util::RequireArg(cws.size() >= required,
      "Dpf::Gen: cws must hold in_bits*stride+1 elements");
  Generate(CwsWriter{cws.data(), stride, tr}, s0s.data(), a, b_buf);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Writer>
__host__ __device__ inline void Dpf<in_bits, Group, Prg, In>::Generate(
    Writer writer, const int4 *s0s, In a, int4 b_buf) {
  int4 s0 = s0s[0];
  s0 = util::SetLsb(s0, false);
  int4 s1 = s0s[1];
  s1 = util::SetLsb(s1, false);
  bool t0 = false;
  bool t1 = true;
  b_buf = util::SetLsb(b_buf, false);
 
  for (int i = 0; i < in_bits; ++i) {
    auto [s0l, s0r] = prg_.Gen(s0);
    auto [s1l, s1r] = prg_.Gen(s1);
 
    bool t0l = util::GetLsb(s0l);
    s0l = util::SetLsb(s0l, false);
    bool t0r = util::GetLsb(s0r);
    s0r = util::SetLsb(s0r, false);
    bool t1l = util::GetLsb(s1l);
    s1l = util::SetLsb(s1l, false);
    bool t1r = util::GetLsb(s1r);
    s1r = util::SetLsb(s1r, false);
 
    bool a_bit = (a >> (in_bits - 1 - i)) & 1;
 
    int4 s_cw;
    if (!a_bit) s_cw = util::Xor(s0r, s1r);
    else s_cw = util::Xor(s0l, s1l);
 
    bool tl_cw = t0l ^ t1l ^ a_bit ^ 1;
    bool tr_cw = t0r ^ t1r ^ a_bit;
 
    if (!a_bit) {
      s0 = s0l;
      if (t0) s0 = util::Xor(s0, s_cw);
      s1 = s1l;
      if (t1) s1 = util::Xor(s1, s_cw);
 
      if (t0) t0 = t0l ^ tl_cw;
      else t0 = t0l;
      if (t1) t1 = t1l ^ tl_cw;
      else t1 = t1l;
    } else {
      s0 = s0r;
      if (t0) s0 = util::Xor(s0, s_cw);
      s1 = s1r;
      if (t1) s1 = util::Xor(s1, s_cw);
 
      if (t0) t0 = t0r ^ tr_cw;
      else t0 = t0r;
      if (t1) t1 = t1r ^ tr_cw;
      else t1 = t1r;
    }
 
    s_cw = util::SetLsb(s_cw, tl_cw);
    writer.Write(i, s_cw, tr_cw);
  }
 
  auto v_cw_np1 = Group::From(b_buf) + (-Group::From(s0)) + Group::From(s1);
  if (t1) v_cw_np1 = -v_cw_np1;
  writer.Finish(v_cw_np1.Into());
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline int4 Dpf<in_bits, Group, Prg, In>::Eval(bool b, int4 s0,
    cuda::std::span<const int4> cws, In tr, In x, size_t stride) {
  util::RequireArg(stride > 0, "Dpf::Eval: stride must be positive");
  constexpr size_t kMaxElements = SIZE_MAX / sizeof(int4);
  util::RequireArg(stride <= (kMaxElements - 1) / static_cast<size_t>(in_bits),
      "Dpf::Eval: stride exceeds the correction array limit");
  const size_t required = static_cast<size_t>(in_bits) * stride + 1;
  util::RequireArg(cws.size() >= required,
      "Dpf::Eval: cws must hold in_bits*stride+1 elements");
  return EvalCws(b, s0, CwsView{cws.data(), tr, stride}, x);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Cws>
__host__ __device__ inline int4 Dpf<in_bits, Group, Prg, In>::EvalCws(
    bool b, int4 s0, Cws cws, In x) {
  int4 state = Walk(b, s0, cws, x, [](bool, bool) {});
  const int4 terminal = cws[in_bits].s;
  util::RequireArg((terminal.w & 1) == 0,
      "Dpf::Eval: the terminal group value must have a zero LSB");
  return ConvertLeaf(b, state, terminal);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
inline void Dpf<in_bits, Group, Prg, In>::EvalRange(bool b, In first, In last,
    int4 s0, cuda::std::span<const int4> cws, In tr, cuda::std::span<int4> ys,
    int par_depth, size_t stride) {
  if (stride == 0) throw std::invalid_argument("Dpf::EvalRange: stride must be positive");
  constexpr size_t kMaxElements = SIZE_MAX / sizeof(int4);
  if (stride > (kMaxElements - 1) / static_cast<size_t>(in_bits))
    throw std::invalid_argument("Dpf::EvalRange: stride exceeds the correction array limit");
  const size_t required = static_cast<size_t>(in_bits) * stride + 1;
  if (cws.size() < required || cws.data() == nullptr)
    throw std::invalid_argument("Dpf::EvalRange: cws must hold in_bits*stride+1 elements");
  EvalRangeCws(b, first, last, s0, CwsView{cws.data(), tr, stride},
      ys, par_depth);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Cws>
inline void Dpf<in_bits, Group, Prg, In>::EvalRangeCws(bool b, In first, In last,
    int4 s0, Cws cws, cuda::std::span<int4> ys, int par_depth) {
  assert(par_depth >= 0);
  if (first > last)
    throw std::invalid_argument("Dpf::EvalRange: first must not exceed last");
  if constexpr (in_bits < sizeof(In) * 8) {
    if ((last >> in_bits) != 0)
      throw std::invalid_argument("Dpf::EvalRange: last must stay inside the input domain");
  }
  const __uint128_t difference = static_cast<__uint128_t>(last) - first;
  if (difference >= std::numeric_limits<size_t>::max() / sizeof(int4))
    throw std::invalid_argument("Dpf::EvalRange: the interval is too large for size_t");
  const size_t count = static_cast<size_t>(difference) + 1;
  if (ys.size() < count || ys.data() == nullptr)
    throw std::invalid_argument("Dpf::EvalRange: ys must hold last-first+1 entries");
  util::RequireArg((cws[in_bits].s.w & 1) == 0,
      "Dpf::EvalRange: the terminal group value must have a zero LSB");
  const size_t partitions = detail::ResolveParPartitions(par_depth, count);
  detail::ParallelFor(par_depth, partitions, [&](size_t partition) {
    const size_t remainder = count % partitions;
    const size_t begin = partition * (count / partitions) +
        (partition < remainder ? partition : remainder);
    const size_t length = count / partitions + (partition < remainder);
    const In partition_first = first + static_cast<In>(begin);
    const In partition_last = partition_first + static_cast<In>(length - 1);
    EvalRangeSerial(b, partition_first, partition_last, s0, cws, ys.data() + begin);
  });
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Cws>
inline void Dpf<in_bits, Group, Prg, In>::EvalRangeSerial(bool b, In first, In last,
    int4 s0, Cws cws, int4 ys[]) {
  constexpr In domain_last = [] {
    if constexpr (in_bits == sizeof(In) * 8) return static_cast<In>(~In{0});
    else return static_cast<In>((In{1} << in_bits) - 1);
  }();
  struct Frame {
    int4 state;
    In low;
    In high;
    int level;
  };
  cuda::std::array<Frame, in_bits + 1> stack;
  size_t pending = 1;
  stack[0] = {util::SetLsb(s0, b), In{0}, domain_last, 0};
  while (pending != 0) {
    const Frame node = stack[--pending];
    if (node.level == in_bits) {
      const int4 terminal = cws[in_bits].s;
      ys[static_cast<size_t>(node.low - first)] = ConvertLeaf(b, node.state, terminal);
      continue;
    }
    const In middle = node.low + ((node.high - node.low) >> 1);
    int4 left;
    int4 right;
    ExpandNode(prg_, node.state, cws[node.level].s, cws[node.level].tr, left, right);
    if (last > middle) {
      assert(pending < stack.size());
      stack[pending++] = {right, static_cast<In>(middle + 1), node.high, node.level + 1};
    }
    if (first <= middle) {
      assert(pending < stack.size());
      stack[pending++] = {left, node.low, middle, node.level + 1};
    }
  }
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
inline void Dpf<in_bits, Group, Prg, In>::EvalAll(bool b, int4 s0,
    cuda::std::span<const int4> cws, In tr, int4 ys[], int par_depth, size_t stride) {
  if (stride == 0) throw std::invalid_argument("Dpf::EvalAll: stride must be positive");
  constexpr size_t kMaxElements = SIZE_MAX / sizeof(int4);
  if (stride > (kMaxElements - 1) / static_cast<size_t>(in_bits))
    throw std::invalid_argument("Dpf::EvalAll: stride exceeds the correction array limit");
  const size_t required = static_cast<size_t>(in_bits) * stride + 1;
  if (cws.size() < required)
    throw std::invalid_argument("Dpf::EvalAll: cws must hold in_bits*stride+1 elements");
  EvalAllCws(b, s0, CwsView{cws.data(), tr, stride}, ys, par_depth);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Cws>
inline void Dpf<in_bits, Group, Prg, In>::EvalAllCws(
    bool b, int4 s0, Cws cws, int4 ys[], int par_depth) {
  assert(par_depth >= 0);
  static_assert(detail::FitsFullDomain<int4, in_bits>(),
      "full-domain output count and byte size must fit in size_t");
 
  // The terminal LSB is a caller precondition. Check it on the calling thread,
  // so a Debug violation aborts there instead of inside an OpenMP worker;
  // Release builds do not check it.
  util::RequireArg((cws[in_bits].s.w & 1) == 0,
      "Dpf::EvalAll: the terminal group value must have a zero LSB");
 
  if constexpr (in_bits < kBatchPrefix + kBatchTail) {
    // Domains too small to fill an 8-wide frontier evaluate per input; the
    // point walk is cheaper than the batched machinery at this size.
    // The Walk call is host-only (it runs inside EvalAllCws's host ParallelFor),
    // so nvcc's #20013-D on the host lambda passed to the host/device Walk
    // template below is a known false alarm for this instantiation.
    constexpr size_t n = size_t{1} << in_bits;
    const int4 terminal = cws[in_bits].s;
    fss::detail::ParallelFor(par_depth, n, [&](size_t x) {
      const int4 state = Walk(b, s0, cws, static_cast<In>(x), [](bool, bool) {});
      ys[x] = ConvertLeaf(b, state, terminal);
    });
    return;
  }
 
  int4 st = s0;
  st = util::SetLsb(st, b);
 
  int par_depth_ = fss::detail::ResolveParDepth(par_depth, in_bits);
 
  int split_depth = par_depth_;
  if (split_depth > in_bits - kBatchPrefix - kBatchTail) split_depth = in_bits - kBatchPrefix - kBatchTail;
  if (split_depth > 0) {
    // Split at a frontier of 2 ** split_depth subtree roots instead of per-node tasks,
    // so the sequential batched walk stays the fast default path.
    size_t num_tasks = size_t{1} << split_depth;
    std::vector<int4> frontier = BuildParallelFrontier(st, cws, split_depth);
    size_t span = (size_t{1} << in_bits) / num_tasks;
    fss::detail::ParallelFor(par_depth_, num_tasks, [&](size_t task) {
      EvalBatched8(b, frontier[task], cws, ys, task * span, split_depth);
    });
  } else {
    EvalBatched8(b, st, cws, ys, 0, 0);
  }
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline void Dpf<in_bits, Group, Prg, In>::ExpandNode(
    Prg &prg, int4 state, int4 correction, bool right_control, int4 &left, int4 &right) {
  int control = static_cast<int>(util::GetLsb(state));
  auto [raw_left, raw_right] = prg.Gen(util::SetLsb(state, false));
  int left_control = static_cast<int>(util::GetLsb(correction));
  int4 seed_correction = util::SetLsb(correction, false);
  left = CorrectChild(raw_left, seed_correction, -control, control & left_control);
  right = CorrectChild(raw_right, seed_correction, -control, control & static_cast<int>(right_control));
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline int4 Dpf<in_bits, Group, Prg, In>::PathStep(
    Prg &prg, int4 state, int4 correction, bool right_control, bool take_right) {
  int control = static_cast<int>(util::GetLsb(state));
  auto [raw_child] = detail::GenHalf<2>(prg, util::SetLsb(state, false), take_right);
  int child_control = static_cast<int>(take_right ? right_control : util::GetLsb(correction));
  return CorrectChild(raw_child, util::SetLsb(correction, false),
      -control, control & child_control);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Correction>
__host__ __device__ inline int4 Dpf<in_bits, Group, Prg, In>::ConvertLeaf(
    bool party, int4 state, Correction &correction) {
  auto value = Group::From(util::SetLsb(state, false));
  if (util::GetLsb(state)) {
    if constexpr (std::is_same_v<std::remove_cv_t<Correction>, int4>) value = value + Group::From(correction);
    else value = value + correction;
  }
  if (party) value = -value;
  return value.Into();
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Cws, typename Visit>
__host__ __device__ inline int4 Dpf<in_bits, Group, Prg, In>::Walk(
    bool b, int4 s0, Cws cws, In x, Visit &&visit) {
  int4 state = util::SetLsb(s0, b);
  for (int i = 0; i < in_bits; ++i) {
    bool x_bit = (x >> (in_bits - 1 - i)) & 1;
    visit(util::GetLsb(state), x_bit);
    state = PathStep(prg_, state, cws[i].s, cws[i].tr, x_bit);
  }
  return state;
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Cws>
inline std::vector<int4> Dpf<in_bits, Group, Prg, In>::BuildParallelFrontier(
    int4 st, Cws cws, int split_depth) {
  std::vector<int4> frontier(size_t{1} << split_depth);
  frontier[0] = st;
  for (int lv = 0; lv < split_depth; ++lv) ExpandFrontier(frontier.data(), size_t{1} << lv, cws, lv);
  return frontier;
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Cws>
inline void Dpf<in_bits, Group, Prg, In>::ExpandFrontier(int4 frontier[], size_t w, Cws cws, int lv) {
  int4 correction = cws[lv].s;
  bool right_control = cws[lv].tr;
  for (size_t j = w; j-- > 0;) {
    int4 left, right;
    ExpandNode(prg_, frontier[j], correction, right_control, left, right);
    frontier[2 * j] = left;
    frontier[2 * j + 1] = right;
  }
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline int4 Dpf<in_bits, Group, Prg, In>::CorrectChild(int4 child, int4 s_cw, int m, int t_cw) {
  int4 res = {child.x ^ (s_cw.x & m), child.y ^ (s_cw.y & m), child.z ^ (s_cw.z & m), child.w ^ (s_cw.w & m)};
  res.w = (res.w & ~1) | ((child.w ^ t_cw) & 1);
  return res;
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Cws>
inline void Dpf<in_bits, Group, Prg, In>::EvalBatched8(
    bool b, int4 st, Cws cws, int4 ys[], size_t l, int i) {
  int4 frontier[8];
  frontier[0] = st;
  for (int lv = i; lv < i + kBatchPrefix; ++lv) ExpandFrontier(frontier, size_t{1} << (lv - i), cws, lv);
  EvalBatched8Group(b, frontier, cws, ys, l, i + kBatchPrefix);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <typename Cws>
inline void Dpf<in_bits, Group, Prg, In>::EvalBatched8Group(
    bool b, const int4 (&frontier)[8], Cws cws, int4 ys[], size_t l, int i) {
  if (i == in_bits - kBatchTail) {
    // Expand the group to its 8 << kBatchTail leaves in local memory, then convert.
    int4 bufs[2][8 << kBatchTail];
    int4 *cur = bufs[0];
    int4 *next = bufs[1];
    for (int j = 0; j < 8; ++j) cur[j] = frontier[j];
    int w = 8;
    for (int lv = i; lv < in_bits; ++lv) {
      int4 s_cw = util::SetLsb(cws[lv].s, false);
      int tl_cw = util::GetLsb(cws[lv].s);
      int tr_cw = cws[lv].tr;
      for (int g = 0; g < w / 8; ++g) {
        cuda::std::array<int4, 8> seeds;
        int ts[8];
        for (int j = 0; j < 8; ++j) {
          ts[j] = static_cast<int>(util::GetLsb(cur[8 * g + j]));
          seeds[j] = util::SetLsb(cur[8 * g + j], false);
        }
        auto out = detail::GenBatch<2>(prg_, seeds);
        for (int j = 0; j < 8; ++j) {
          int m = -ts[j];
          next[16 * g + 2 * j] = CorrectChild(out[j][0], s_cw, m, ts[j] & tl_cw);
          next[16 * g + 2 * j + 1] = CorrectChild(out[j][1], s_cw, m, ts[j] & tr_cw);
        }
      }
      int4 *swap = cur;
      cur = next;
      next = swap;
      w *= 2;
    }
 
    int4 v_cw_np1 = cws[in_bits].s;
    // The shared entry validated the terminal before parallel dispatch.
    assert((v_cw_np1.w & 1) == 0);
    Group v_cw = Group::From(v_cw_np1);
    for (int j = 0; j < w; ++j) ys[l + j] = ConvertLeaf(b, cur[j], v_cw);
    return;
  }
 
  int4 s_cw = util::SetLsb(cws[i].s, false);
  int tl_cw = util::GetLsb(cws[i].s);
  int tr_cw = cws[i].tr;
 
  cuda::std::array<int4, 8> seeds;
  int ts[8];
  for (int j = 0; j < 8; ++j) {
    ts[j] = static_cast<int>(util::GetLsb(frontier[j]));
    seeds[j] = util::SetLsb(frontier[j], false);
  }
  auto out = detail::GenBatch<2>(prg_, seeds);
  // The 16 children interleave in leaf order: node j's children land at 2j, 2j+1 of the
  // next level. De-interleave them into two consecutive groups of 8.
  int4 low[8];
  int4 high[8];
  for (int j = 0; j < 4; ++j) {
    int m = -ts[j];
    low[2 * j] = CorrectChild(out[j][0], s_cw, m, ts[j] & tl_cw);
    low[2 * j + 1] = CorrectChild(out[j][1], s_cw, m, ts[j] & tr_cw);
    int m4 = -ts[j + 4];
    high[2 * j] = CorrectChild(out[j + 4][0], s_cw, m4, ts[j + 4] & tl_cw);
    high[2 * j + 1] = CorrectChild(out[j + 4][1], s_cw, m4, ts[j + 4] & tr_cw);
  }
 
  size_t half = size_t{8} << (in_bits - i - 1);
  EvalBatched8Group(b, low, cws, ys, l, i + 1);
  EvalBatched8Group(b, high, cws, ys, l + half, i + 1);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <int frontier_depth, int threads_per_block>
struct Dpf<in_bits, Group, Prg, 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 constexpr int kFrontierDepth = frontier_depth < in_bits ? frontier_depth : in_bits;
  static_assert(in_bits - kFrontierDepth <= 8,
      "gpu frontier leaves more than eight tail levels; specify a larger frontier_depth");
  static constexpr int kThreadsPerBlockLog2 = [] {
    int threads_per_block_log2 = 0;
    for (int thread_count = threads_per_block; thread_count > 1; thread_count >>= 1) ++threads_per_block_log2;
    return threads_per_block_log2 < kFrontierDepth ? threads_per_block_log2 : 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, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <int frontier_depth, int threads_per_block>
__host__ inline cudaError_t Dpf<in_bits, Group, Prg, In>::EvalAllGpuBatchCw(
    bool b, cuda::std::span<const int4> s0s, cuda::std::span<const Cw> cws,
    cuda::std::span<int4> ys, cudaStream_t stream) const {
#if defined(__CUDACC__)
  using LaunchShape = GpuLaunchShape<frontier_depth, threads_per_block>;
  using Scheme = Dpf<in_bits, Group, Prg, In>;
  constexpr size_t kOutputsPerKey = size_t{1} << in_bits;
  constexpr size_t kCwPerKey = in_bits + 1;
  constexpr size_t kMaxAddressableSize = std::numeric_limits<size_t>::max();
  size_t key_count = s0s.size();
  if (key_count == 0) return cudaSuccess;
  if (key_count > kMaxAddressableSize / kCwPerKey / sizeof(Cw) || key_count > kMaxAddressableSize / kOutputsPerKey / sizeof(int4) ||
      cws.size() < key_count * kCwPerKey || 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::DpfEvalAllBatchKernel<Scheme, LaunchShape::kFrontierDepth, LaunchShape::kBlockRootDepth>
      <<<key_count * LaunchShape::kBlocksPerKey, LaunchShape::kThreadsPerBlock, 0, stream>>>(
          b, s0s.data(), cws.data(), ys.data(), prg_);
  return cudaGetLastError();
#else
  static_assert(sizeof(Prg) == 0, "gpu launch requires a CUDA compiler");
  return cudaErrorNotSupported;
#endif
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <int frontier_depth, int threads_per_block>
__host__ inline cudaError_t Dpf<in_bits, Group, Prg, In>::EvalRangeGpu(
    bool b, In first, In last, int4 s0, cuda::std::span<const int4> cws, In tr,
    cuda::std::span<int4> ys, cudaStream_t stream, size_t stride) const {
#if defined(__CUDACC__)
  using Scheme = Dpf<in_bits, Group, Prg, In>;
  static_assert(frontier_depth >= 0, "gpu frontier depth must be nonnegative");
  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]");
  constexpr size_t kMaxElements = std::numeric_limits<size_t>::max() / sizeof(int4);
  if (stride == 0 || (stride > (kMaxElements - 1) / static_cast<size_t>(in_bits)) ||
      cws.size() < static_cast<size_t>(in_bits) * stride + 1 ||
      cws.data() == nullptr || ys.data() == nullptr)
    return cudaErrorInvalidValue;
  return detail::SplitRangeChunks<in_bits, frontier_depth, threads_per_block>(
      first, last, ys.size(),
      [&](In chunk_first, int depth, size_t offset, int frontier, int thread_bits, int root,
          size_t blocks) {
        gpu::detail::DpfSplitEvalRangeKernel<Scheme, threads_per_block>
            <<<blocks, 1 << thread_bits, 0, stream>>>(b, s0, cws.data(), tr, stride,
                ys.data() + offset, chunk_first, depth, frontier, root, prg_);
        return cudaGetLastError();
      });
#else
  static_assert(sizeof(Prg) == 0, "gpu launch requires a CUDA compiler");
  return cudaErrorNotSupported;
#endif
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <int frontier_depth, int threads_per_block>
__host__ inline cudaError_t Dpf<in_bits, Group, Prg, In>::EvalAllGpu(
    bool b, int4 s0, cuda::std::span<const int4> cws, In tr,
    cuda::std::span<int4> ys, cudaStream_t stream, size_t stride) const
  requires(in_bits < sizeof(size_t) * 8) {
#if defined(__CUDACC__)
  using LaunchShape = GpuLaunchShape<frontier_depth, threads_per_block>;
  using Scheme = Dpf<in_bits, Group, Prg, In>;
  constexpr size_t kOutputsPerKey = size_t{1} << in_bits;
  constexpr size_t kMaxElements = std::numeric_limits<size_t>::max() / sizeof(int4);
  if (stride == 0 ||
      (stride > (kMaxElements - 1) / static_cast<size_t>(in_bits)))
    return cudaErrorInvalidValue;
  const size_t required = static_cast<size_t>(in_bits) * 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::DpfSplitEvalAllKernel<Scheme, LaunchShape::kFrontierDepth, LaunchShape::kBlockRootDepth>
      <<<LaunchShape::kBlocksPerKey, LaunchShape::kThreadsPerBlock, 0, stream>>>(
          b, s0, cws.data(), tr, stride, ys.data(), prg_);
  return cudaGetLastError();
#else
  static_assert(sizeof(Prg) == 0, "gpu launch requires a CUDA compiler");
  return cudaErrorNotSupported;
#endif
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
template <int frontier_depth, int threads_per_block>
__host__ inline cudaError_t Dpf<in_bits, Group, Prg, In>::EvalAllGpuBatch(
    bool b, cuda::std::span<const int4> s0s, cuda::std::span<const int4> cws,
    cuda::std::span<const In> tr, cuda::std::span<int4> ys,
    cudaStream_t stream, size_t stride) const
  requires(in_bits < sizeof(size_t) * 8) {
#if defined(__CUDACC__)
  using LaunchShape = GpuLaunchShape<frontier_depth, threads_per_block>;
  using Scheme = Dpf<in_bits, Group, Prg, In>;
  constexpr size_t kOutputsPerKey = size_t{1} << in_bits;
  constexpr size_t kMaxElements = std::numeric_limits<size_t>::max() / sizeof(int4);
  constexpr size_t kMaxAddressableSize = std::numeric_limits<size_t>::max();
  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>(in_bits))
    return cudaErrorInvalidValue;
  const size_t required = static_cast<size_t>(in_bits) * resolved_stride + key_count;
  if (cws.size() < required || tr.size() < key_count ||
      ys.size() < key_count * kOutputsPerKey ||
      s0s.data() == nullptr || cws.data() == nullptr || tr.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::DpfSplitEvalAllBatchKernel<Scheme, LaunchShape::kFrontierDepth, LaunchShape::kBlockRootDepth>
      <<<key_count * LaunchShape::kBlocksPerKey, LaunchShape::kThreadsPerBlock, 0, stream>>>(
          b, s0s.data(), cws.data(), tr.data(), resolved_stride, ys.data(), prg_);
  return cudaGetLastError();
#else
  static_assert(sizeof(Prg) == 0, "gpu launch requires a CUDA compiler");
  return cudaErrorNotSupported;
#endif
}
 
#if FSS_ENABLE_V1_COMPAT
 
// Deprecated v1.3.0 aggregate-Cw entries. Every wrapper copies between the
// caller's Cw layout and a compile-time-sized split-component scratch buffer,
// then runs the current public method at stride 1. The scratch copies are the
// compatibility cost and never enter the core algorithms.
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline void Dpf<in_bits, Group, Prg, In>::Gen(
    Cw cws[], const int4 s0s[2], In a, int4 b_buf) {
  int4 seed[in_bits + 1];
  In tr;
  Gen(cuda::std::span<int4>(seed, in_bits + 1), tr,
      cuda::std::span<const int4, 2>(s0s, 2), a, b_buf);
  for (int i = 0; i < in_bits; ++i) cws[i] = {seed[i], ((tr >> i) & 1) != 0};
  // The terminal is a plain group value; its unused tr bit is false.
  cws[in_bits] = {seed[in_bits], false};
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline int4 Dpf<in_bits, Group, Prg, In>::Eval(
    bool b, int4 s0, const Cw cws[], In x) {
  int4 seed[in_bits + 1];
  In tr = 0;
  for (int i = 0; i <= in_bits; ++i) seed[i] = cws[i].s;
  for (int i = 0; i < in_bits; ++i) tr |= static_cast<In>(cws[i].tr) << i;
  return Eval(b, s0, cuda::std::span<const int4>(seed, in_bits + 1), tr, x);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
inline void Dpf<in_bits, Group, Prg, In>::EvalAll(
    bool b, int4 s0, const Cw cws[], int4 ys[]) {
  int4 seed[in_bits + 1];
  In tr = 0;
  for (int i = 0; i <= in_bits; ++i) seed[i] = cws[i].s;
  for (int i = 0; i < in_bits; ++i) tr |= static_cast<In>(cws[i].tr) << i;
  EvalAll(b, s0, cuda::std::span<const int4>(seed, in_bits + 1), tr, ys);
}
 
// Deprecated v1.3.0 aggregate-Key entries. Same conversion as the Cw entries
// above; the Key entries additionally spread the shared correction words over
// both keys and take the party seeds from s0s. The bodies spell the current
// public Key type, which stays undeprecated.
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline void Dpf<in_bits, Group, Prg, In>::Gen(
    Key &k0, Key &k1, cuda::std::span<const int4, 2> s0s, In a, int4 b_buf) {
  int4 seed[in_bits + 1];
  In tr;
  Gen(cuda::std::span<int4>(seed, in_bits + 1), tr, s0s, a, b_buf);
  for (int i = 0; i < in_bits; ++i) {
    k0.cws[i] = {seed[i], ((tr >> i) & 1) != 0};
    k1.cws[i] = k0.cws[i];
  }
  // The terminal is a plain group value; its unused tr bit is false.
  k0.cws[in_bits] = {seed[in_bits], false};
  k1.cws[in_bits] = k0.cws[in_bits];
  k0.s0 = s0s[0];
  k1.s0 = s0s[1];
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
__host__ __device__ inline int4 Dpf<in_bits, Group, Prg, In>::Eval(
    bool b, const Key &key, In x) {
  int4 seed[in_bits + 1];
  In tr = 0;
  for (int i = 0; i <= in_bits; ++i) seed[i] = key.cws[i].s;
  for (int i = 0; i < in_bits; ++i) tr |= static_cast<In>(key.cws[i].tr) << i;
  return Eval(b, key.s0, cuda::std::span<const int4>(seed, in_bits + 1), tr, x);
}
 
template <int in_bits, typename Group, typename Prg, typename In>
  requires(ValidForInputBits<In, in_bits> &&
      Groupable<Group> && Prgable<Prg, 2>)
inline void Dpf<in_bits, Group, Prg, In>::EvalAll(
    bool b, const Key &key, cuda::std::span<int4> ys) {
  int4 seed[in_bits + 1];
  In tr = 0;
  for (int i = 0; i <= in_bits; ++i) seed[i] = key.cws[i].s;
  for (int i = 0; i < in_bits; ++i) tr |= static_cast<In>(key.cws[i].tr) << i;
  EvalAll(b, key.s0, cuda::std::span<const int4>(seed, in_bits + 1), tr, ys.data());
}
 
#endif  // FSS_ENABLE_V1_COMPAT
 
}  // namespace fss