/// @file dpf/eval_sequence.hpp
/// @brief Evaluate a sorted list of DPF inputs.
/// @details The range is nondecreasing; an unsorted range throws
///          `std::runtime_error`. `return_output_only_tag_` stores one share
///          per point. `return_entire_node_tag_` stores whole leaves and is
///          the default. A `sequence_recipe` repeats the list, and a sequence
///          memoizer bound to that recipe object resumes the traversal.
///          @snippet evaluation/eval_sequence.cpp eval-sequence
/// @author Ryan Henry <ryan.henry@ucalgary.ca>
/// @author Christopher Jiang <christopher.jiang@ucalgary.ca>
/// @copyright Copyright (c) 2019-2024 Ryan Henry and [others](@ref authors)
/// @license Released under a GNU General Public v2.0 (GPLv2) license;
///          see [LICENSE.md](@ref license) for details.

#ifndef LIBDPF_INCLUDE_DPF_EVAL_SEQUENCE_HPP__
#define LIBDPF_INCLUDE_DPF_EVAL_SEQUENCE_HPP__

#include <portable-snippets/builtin/builtin.h>
#include "hedley/hedley.h"

#include <cstddef>
#include <cstring>
#include <type_traits>
#include <utility>
#include <tuple>
#include <algorithm>
#include <array>
#include <iterator>
#include <stdexcept>
#include <limits>
#include <vector>

#include "dpf/dpf_key.hpp"
#include "dpf/eval_common.hpp"
#include "dpf/eval_target.hpp"
#include "dpf/eval_point.hpp"
#include "dpf/eval_interval.hpp"
#include "dpf/path_memoizer.hpp"
#include "dpf/sequence_memoizer.hpp"
#include "dpf/sequence_utils.hpp"
#include "dpf/subsequence_iterable.hpp"
#include "dpf/subinterval_iterable.hpp"
#include "dpf/verifiable.hpp"

namespace dpf
{

template <typename KeyT, typename ForwardIterator>
void prove_sequence(const KeyT & key, ForwardIterator begin, ForwardIterator end,
    prove_ref pr);

// Contiguous runs at least this long use the interval tree (one expand per
// node). Shorter runs stay on the path memoizer. Isolated points at least
// this many use the breadth-first block walk, which beats a fresh path per
// point once the list is wide.
inline constexpr std::size_t sequence_interval_run = 24;
// Breadth-first shares prefixes on sparse lists, but its per-level block-split
// bookkeeping loses to a path memoizer once single-child expands cut the
// path PRG cost in half. Keep the explicit `eval_sequence_breadth_first`
// entry point; do not auto-select it from `eval_sequence`.
inline constexpr std::size_t sequence_breadth_points = std::numeric_limits<std::size_t>::max();

template <std::size_t I,
          typename DpfKey,
          typename ForwardIterator,
          typename OutputBuffer>
inline auto eval_sequence_breadth_first(const DpfKey & dpf, ForwardIterator begin,
    ForwardIterator end, OutputBuffer && outbuf);

template <bool Entire,
          std::size_t ...Is,
          typename DpfKey,
          typename InputT,
          typename OutputBuffers,
          std::size_t ...IIs>
void scatter_interval_run(const DpfKey & dpf, InputT from, InputT to,
    OutputBuffers & outbufs, std::size_t index0, std::index_sequence<IIs...>)
{
    // std::make_tuple, not utils::make_tuple: one output must stay a tuple
    // so each selected buffer is addressable by index.
    auto ibufs = std::make_tuple(
        make_output_buffer_for_interval<DpfKey, Is>(from, to)...);
    (void)eval_interval<Is...>(dpf, from, to, ibufs);

    constexpr auto opl = DpfKey::outputs_per_leaf;
    constexpr auto lg = DpfKey::lg_outputs_per_leaf;
    constexpr auto mod_pow_2 = utils::mod_pow_2<InputT>{};
    const std::size_t preclip = mod_pow_2(dpf.offset_x(from), lg);

    constexpr auto to_int = utils::to_integral_type<InputT>{};
    auto span = to_int(to) - to_int(from);
    constexpr auto bits = utils::bitlength_of_v<InputT>;
    if constexpr (bits < utils::bitlength_of_v<decltype(span)>)
        span &= (decltype(span){1} << bits) - 1;
    const std::size_t n = static_cast<std::size_t>(span) + 1;
    constexpr std::size_t leaf_mask = (std::size_t{1} << lg) - 1;

    auto one = [&](auto which)
    {
        constexpr std::size_t k = decltype(which)::value;
        auto & src = utils::get<k>(ibufs);
        auto & dst = utils::get<k>(outbufs);
        auto put = [](auto & slot, const auto & val)
        {
            assign_share_slot(slot, val);
        };
        if constexpr (Entire)
        {
            for (std::size_t j = 0; j < n; ++j)
            {
                const std::size_t elem_i = preclip + j;
                const std::size_t leaf = elem_i - (elem_i & leaf_mask);
                for (std::size_t lane = 0; lane < opl; ++lane)
                    put(dst[(index0 + j) * opl + lane], src[leaf + lane]);
            }
        }
        else
        {
            for (std::size_t j = 0; j < n; ++j)
                put(dst[index0 + j], src[preclip + j]);
        }
    };
    (one(std::integral_constant<std::size_t, IIs>{}), ...);
}

namespace internal
{

/// @brief Leaf interior nodes for a list, eight paths at a time.
/// @details One scalar spine up to the shallowest divergence, then
///          `tree::traverse8` across the whole suffix of each level. `walk[i]`
///          is the offset input with the signed MSB already flipped, matching
///          `eval_point`.
template <typename DpfKey>
std::vector<typename DpfKey::interior_node> sequence_wide_leaves(
    const DpfKey & dpf, const typename DpfKey::input_type * walk, std::size_t n)
{
    using node = typename DpfKey::interior_node;
    std::vector<node> leaves(n);
    if (n == 0)
        return leaves;
    constexpr auto clz = utils::countl_zero_symmetric_difference<
        typename DpfKey::input_type>{};
    std::size_t common = dpf.depth + 1;
    for (std::size_t i = 1; i < n; ++i)
        common = std::min(common, clz(walk[0], walk[i]) + 1);

    std::array<node, DpfKey::depth + 1> path{};
    path[0] = dpf.root();
    auto mask = dpf.msb_mask;
    for (std::size_t level = 1; level < common && level <= dpf.depth;
         ++level, mask >>= 1)
    {
        const bool bit = (mask & walk[0]) != 0;
        path[level] = DpfKey::traverse_interior(path[level - 1],
            dpf.correction_word(level - 1, bit), bit,
            DpfKey::tree::is_last_level(level - 1, dpf.depth));
    }
    if (common > dpf.depth)
    {
        std::fill(leaves.begin(), leaves.end(), path[dpf.depth]);
        return leaves;
    }

    std::vector<node> cur(n, path[common - 1]);
    std::vector<node> next(n);
    mask = dpf.msb_mask >> (common - 1);
    for (std::size_t level = common; level <= dpf.depth; ++level, mask >>= 1)
    {
        const node cw0 = dpf.correction_word(level - 1, false);
        const node cw1 = dpf.correction_word(level - 1, true);
        const bool last = DpfKey::tree::is_last_level(level - 1, dpf.depth);
        std::size_t i = 0;
        for (; i + 8 <= n; i += 8)
        {
            node parents[8];
            node cws[8];
            bool dirs[8];
            for (int k = 0; k < 8; ++k)
            {
                dirs[k] = (mask & walk[i + static_cast<std::size_t>(k)]) != 0;
                cws[k] = dirs[k] ? cw1 : cw0;
                parents[k] = cur[i + static_cast<std::size_t>(k)];
            }
            node outs[8];
            DpfKey::tree::traverse8(parents, cws, dirs, outs);
            for (int k = 0; k < 8; ++k)
                next[i + static_cast<std::size_t>(k)] = outs[k];
        }
        for (; i < n; ++i)
        {
            const bool bit = (mask & walk[i]) != 0;
            next[i] = DpfKey::traverse_interior(cur[i], bit ? cw1 : cw0, bit, last);
        }
        cur.swap(next);
    }
    return cur;
}

template <std::size_t ...Is,
          typename DpfKey,
          typename ForwardIterator,
          typename OutputBuffers,
          std::size_t ...IIs>
auto eval_sequence_entire_node(const DpfKey & dpf, ForwardIterator begin, ForwardIterator end,
    OutputBuffers && outbufs, std::index_sequence<IIs...>)
{
    static constexpr std::size_t outputs_per_leaf = DpfKey::outputs_per_leaf;
    using input_type = typename DpfKey::input_type;
    constexpr bool any_packed =
        (utils::is_packed_subbyte_v<typename DpfKey::concrete_output_type<Is>> || ...);
    constexpr bool any_bit =
        (std::is_same_v<typename DpfKey::concrete_output_type<Is>, dpf::bit> || ...);
    constexpr bool integral_in = std::is_integral_v<input_type>;
    // Bit and packed leaves are proxy buffers. The interval scatter writes
    // ordinary word slots; those outputs stay on the path memoizer.
    constexpr bool interval_ok = integral_in && !any_packed && !any_bit;

    auto write_point = [&](auto & path, std::size_t i, const input_type & x)
    {
        if constexpr (any_packed)
        {
            auto nodes = std::make_tuple(dpf::eval_point<Is>(dpf, x, path).node...);
            (store_leaf_bytes(utils::get<IIs>(outbufs), i, std::get<IIs>(nodes)), ...);
        }
        else
        {
            auto temp = std::make_tuple(dpf::eval_point<Is>(dpf, x, path).node...);
            (utils::raw_memcpy(&utils::get<IIs>(outbufs)[i * outputs_per_leaf],
                &utils::get<IIs>(temp),
                sizeof(typename DpfKey::concrete_output_type<Is>) * outputs_per_leaf), ...);
        }
    };

    if constexpr (interval_ok)
    {
        using iter_cat = typename std::iterator_traits<ForwardIterator>::iterator_category;
        bool classify = true;
        if constexpr (std::is_base_of_v<std::random_access_iterator_tag, iter_cat>)
        {
            // Below this length the interval-run cover does not apply, so skip
            // the classification pass.
            classify = static_cast<std::size_t>(std::distance(begin, end))
                >= sequence_interval_run;
        }
        if (classify)
        {
        bool sorted = true;
        std::size_t n = 0;
        std::size_t longest = 0;
        std::size_t cur = 0;
        bool dup = false;
        input_type prev{};
        for (auto it = begin; it != end; ++it, ++n)
        {
            if (n == 0)
            {
                prev = *it;
                longest = 1;
                cur = 1;
                continue;
            }
            const input_type x = *it;
            if (x < prev)
                sorted = false;
            else if (x == prev)
                dup = true;
            if (x == static_cast<input_type>(prev + input_type{1}))
            {
                ++cur;
                if (cur > longest)
                    longest = cur;
            }
            else
                cur = 1;
            prev = x;
        }
        if (sorted && longest >= sequence_interval_run)
        {
            std::size_t i = 0;
            for (auto it = begin; it != end; )
            {
                const input_type run_from = *it;
                input_type run_to = run_from;
                auto run_end = it;
                ++run_end;
                while (run_end != end
                    && *run_end == static_cast<input_type>(run_to + input_type{1}))
                {
                    run_to = *run_end;
                    ++run_end;
                }
                const auto span_n = static_cast<std::size_t>(std::distance(it, run_end));
                if (span_n >= sequence_interval_run)
                {
                    scatter_interval_run<true, Is...>(dpf, run_from, run_to, outbufs, i,
                        std::index_sequence<IIs...>{});
                    i += span_n;
                }
                else
                {
                    auto path = make_basic_path_memoizer(dpf);
                    for (auto p = it; p != run_end; ++p, ++i)
                        write_point(path, i, *p);
                }
                it = run_end;
            }
            return utils::make_tuple(
                dpf::subsequence_iterable<DpfKey, decltype(std::begin(utils::get<IIs>(outbufs))), ForwardIterator>(std::begin(utils::get<IIs>(outbufs)), begin, end)...);
        }
        // Breadth-first pays for a block split at every level. On a short
        // tree (8-bit inputs) that overhead loses to the path memoizer.
        // The one-output overload is the only one; keep it out of the
        // multi-output instantiation.
        if constexpr (sizeof...(Is) == 1)
        {
        if (sorted && !dup && n >= sequence_breadth_points
            && DpfKey::depth >= 16)
        {
            eval_sequence_breadth_first<Is...>(dpf, begin, end, utils::get<0>(outbufs));
            return utils::make_tuple(
                dpf::subsequence_iterable<DpfKey, decltype(std::begin(utils::get<IIs>(outbufs))), ForwardIterator>(std::begin(utils::get<IIs>(outbufs)), begin, end)...);
        }
        }
        }
    }

    const std::size_t nseq = static_cast<std::size_t>(std::distance(begin, end));
    if constexpr (sizeof...(Is) == 1 && !DpfKey::is_verifiable
        && prg_has_indep4<typename DpfKey::interior_prg>::value)
    {
        if (nseq >= 16)
        {
            using input_type = typename DpfKey::input_type;
            std::vector<input_type> walk(nseq);
            std::vector<input_type> lane(nseq);
            std::size_t i = 0;
            for (auto it = begin; it != end; ++it, ++i)
            {
                lane[i] = dpf.offset_x(*it);
                walk[i] = lane[i];
                utils::flip_msb_if_signed_integral(walk[i]);
            }
            const auto leaves = sequence_wide_leaves(dpf, walk.data(), nseq);
            constexpr std::size_t ids[] = {Is...};
            constexpr std::size_t I0 = ids[0];
            using output_type = typename DpfKey::concrete_output_type<I0>;
            auto & dst = utils::get<0>(outbufs);
            for (std::size_t p = 0; p < nseq; ++p)
            {
                auto wrapped = make_eval_dpf_output<DpfKey, output_type>(
                    dpf.template traverse_exterior<I0>(leaves[p]), lane[p]);
                utils::raw_memcpy(&dst[p * outputs_per_leaf], &wrapped.node,
                    sizeof(output_type) * outputs_per_leaf);
            }
            return utils::make_tuple(
                dpf::subsequence_iterable<DpfKey, decltype(std::begin(utils::get<IIs>(outbufs))), ForwardIterator>(std::begin(utils::get<IIs>(outbufs)), begin, end)...);
        }
    }

    auto path = make_basic_path_memoizer(dpf);
    std::size_t i = 0;
    for (auto it = begin; it != end; ++it, ++i)
        write_point(path, i, *it);
    return utils::make_tuple(
        dpf::subsequence_iterable<DpfKey, decltype(std::begin(utils::get<IIs>(outbufs))), ForwardIterator>(std::begin(utils::get<IIs>(outbufs)), begin, end)...);
}

template <typename Slot, typename Val>
HEDLEY_ALWAYS_INLINE
void assign_eval_slot(Slot && slot, Val && val)
{
    assign_share_slot(slot, std::forward<Val>(val));
}

template <std::size_t ...Is,
          typename DpfKey,
          typename ForwardIterator,
          typename OutputBuffers,
          std::size_t ...IIs>
auto eval_sequence_output_only(const DpfKey & dpf, ForwardIterator begin, ForwardIterator end, OutputBuffers && outbufs,
    std::index_sequence<IIs...>)
{
    using input_type = typename DpfKey::input_type;
    constexpr bool any_packed =
        (utils::is_packed_subbyte_v<typename DpfKey::concrete_output_type<Is>> || ...);
    constexpr bool any_bit =
        (std::is_same_v<typename DpfKey::concrete_output_type<Is>, dpf::bit> || ...);
    constexpr bool integral_in = std::is_integral_v<input_type>;

    auto write_point = [&](auto & path, std::size_t i, const input_type & x)
    {
        (assign_eval_slot(utils::get<IIs>(outbufs)[i],
             *dpf::eval_point<Is>(dpf, x, path)), ...);
    };

    auto finish = [&](std::size_t i)
    {
        if (i == 0)
        {
            return utils::make_tuple(subinterval_iterable(std::begin(utils::get<IIs>(outbufs)), utils::size(utils::get<IIs>(outbufs)), 0, 0, 0, 0, false)...);
        }
        return utils::make_tuple(subinterval_iterable(std::begin(utils::get<IIs>(outbufs)), utils::size(utils::get<IIs>(outbufs)), 0, i - 1, 0, 0)...);
    };

    if constexpr (integral_in && !any_packed && !any_bit)
    {
        bool sorted = true;
        std::size_t longest = 0;
        std::size_t cur = 0;
        std::size_t n = 0;
        input_type prev{};
        for (auto it = begin; it != end; ++it, ++n)
        {
            if (n == 0)
            {
                prev = *it;
                longest = 1;
                cur = 1;
                continue;
            }
            const input_type x = *it;
            if (x < prev)
                sorted = false;
            if (x == static_cast<input_type>(prev + input_type{1}))
            {
                ++cur;
                if (cur > longest)
                    longest = cur;
            }
            else
                cur = 1;
            prev = x;
        }
        if (sorted && longest >= sequence_interval_run)
        {
            auto path = make_basic_path_memoizer(dpf);
            std::size_t i = 0;
            for (auto it = begin; it != end; )
            {
                const input_type run_from = *it;
                input_type run_to = run_from;
                auto run_end = it;
                ++run_end;
                while (run_end != end
                    && *run_end == static_cast<input_type>(run_to + input_type{1}))
                {
                    run_to = *run_end;
                    ++run_end;
                }
                const auto span_n = static_cast<std::size_t>(std::distance(it, run_end));
                if (span_n >= sequence_interval_run)
                {
                    scatter_interval_run<false, Is...>(dpf, run_from, run_to, outbufs, i,
                        std::index_sequence<IIs...>{});
                    i += span_n;
                }
                else
                {
                    for (auto p = it; p != run_end; ++p, ++i)
                        write_point(path, i, *p);
                }
                it = run_end;
            }
            return finish(i);
        }
    }

    auto path = make_basic_path_memoizer(dpf);
    std::size_t i = 0;
    for (auto it = begin; it != end; ++it, ++i)
        write_point(path, i, *it);
    return finish(i);
}

}  // namespace internal

/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          std::size_t ...Is,
          typename DpfKey,
          typename ForwardIterator,
          typename OutputBuffers,
          typename ReturnType = return_entire_node_tag_,
          std::enable_if_t<looks_like_dpf_key_v<DpfKey> && !is_multilevel_key_v<DpfKey>, bool> = true,
          std::enable_if_t<std::is_base_of_v<return_type_tag_, ReturnType>, bool> = true,
          std::enable_if_t<!std::is_base_of_v<return_type_tag_, OutputBuffers>, bool> = true>
inline auto eval_sequence(const DpfKey & dpf, ForwardIterator begin, ForwardIterator end,
    OutputBuffers && outbufs, ReturnType return_type = ReturnType{})
{
    (void)return_type;
    static_assert(std::is_same_v<ReturnType, return_entire_node_tag_> ||
                    std::is_same_v<ReturnType, return_output_only_tag_>);
    if constexpr(std::is_same_v<ReturnType, return_entire_node_tag_>)
    {
        return internal::eval_sequence_entire_node<I, Is...>(dpf, begin, end, outbufs, std::make_index_sequence<1+sizeof...(Is)>{});
    }
    else
    {
        return internal::eval_sequence_output_only<I, Is...>(dpf, begin, end, outbufs, std::make_index_sequence<1+sizeof...(Is)>{});
    }
}

/// @brief Evaluate a sorted sequence and fold the same VDPF proof as `prove_sequence`.
/// @details Uses interval-run covers (not a path-memo double walk). Values are
///          written by a single subsequent `eval_sequence` without proof.
/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          std::size_t ...Is,
          typename DpfKey,
          typename ForwardIterator,
          typename OutputBuffers,
          typename ReturnType = return_entire_node_tag_,
          std::enable_if_t<looks_like_dpf_key_v<DpfKey> && !is_multilevel_key_v<DpfKey>, bool> = true,
          std::enable_if_t<std::is_base_of_v<return_type_tag_, ReturnType>, bool> = true,
          std::enable_if_t<!std::is_base_of_v<return_type_tag_, OutputBuffers>, bool> = true>
inline auto eval_sequence(const DpfKey & dpf, ForwardIterator begin, ForwardIterator end,
    OutputBuffers && outbufs, prove_ref pr, ReturnType return_type = ReturnType{})
{
    static_assert(DpfKey::is_verifiable,
        "eval_sequence(..., prove(π)): key must carry dpf::verifiable");
    prove_sequence(dpf, begin, end, pr);
    return eval_sequence<I, Is...>(dpf, begin, end,
        std::forward<OutputBuffers>(outbufs), return_type);
}

/// @brief Evaluate a sorted sequence and fold each written output into a sketch.
/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          std::size_t ...Is,
          typename DpfKey,
          typename ForwardIterator,
          typename OutputBuffers,
          typename ReturnType = return_entire_node_tag_,
          std::enable_if_t<looks_like_dpf_key_v<DpfKey> && !is_multilevel_key_v<DpfKey>, bool> = true,
          std::enable_if_t<std::is_base_of_v<return_type_tag_, ReturnType>, bool> = true,
          std::enable_if_t<!std::is_base_of_v<return_type_tag_, OutputBuffers>, bool> = true>
inline auto eval_sequence(const DpfKey & dpf, ForwardIterator begin, ForwardIterator end,
    OutputBuffers && outbufs, sketch_ref & sk, ReturnType return_type = ReturnType{})
{
    static_assert(DpfKey::is_extractable,
        "eval_sequence(..., sketch(σ)): key must carry dpf::extractable");
    auto ret = eval_sequence<I, Is...>(dpf, begin, end, outbufs, return_type);
    if constexpr (sizeof...(Is) == 0)
    {
        for (std::size_t k = 0; k < utils::size(outbufs); ++k)
            sk.absorb(outbufs[k]);
    }
    return ret;
}

/// @brief Evaluate the sorted range `[begin, end)`, allocating a buffer.
/// @tparam I output index
/// @tparam Is is
/// @tparam DpfKey DPF key type
/// @tparam ForwardIterator forward iterator type
/// @tparam ReturnType return type
/// @tparam DpfKey DPF key type
/// @tparam ReturnType return type
/// @param return_type `return_entire_node_tag_{}` or `return_output_only_tag_{}`.
/// @param dpf the DPF key
/// @param begin the iterator to the first query
/// @param end the iterator past the last query
/// @return Pair of buffer (or tuple of buffers) and an iterable in list order.
/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          std::size_t ...Is,
          typename DpfKey,
          typename ForwardIterator,
          typename ReturnType = return_entire_node_tag_,
          std::enable_if_t<looks_like_dpf_key_v<DpfKey> && !is_multilevel_key_v<DpfKey>, bool> = true,
          std::enable_if_t<std::is_base_of_v<return_type_tag_, ReturnType>, bool> = true>
auto eval_sequence(const DpfKey & dpf, ForwardIterator begin, ForwardIterator end,
    ReturnType return_type = ReturnType{})
{
    auto outbufs = utils::make_tuple(
        make_output_buffer_for_subsequence<I>(dpf, begin, end, return_type),
        make_output_buffer_for_subsequence<Is>(dpf, begin, end, return_type)...);

    // moving `outbufs` is allowed as the `outbufs` are `std::vectors`
    //   the underlying data remains on the heap
    //   and thus the data the iterable refers to is still valid
    auto iterable = eval_sequence<I, Is...>(dpf, begin, end, outbufs, return_type);
    return std::make_pair(std::move(outbufs), std::move(iterable));
}

/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          typename DpfKey,
          typename ForwardIterator,
          typename OutputBuffer>
inline auto eval_sequence_breadth_first(const DpfKey & dpf, ForwardIterator begin, ForwardIterator end, OutputBuffer && outbuf)
{
    assert_not_wildcard_output<I>(dpf);

    using dpf_type = DpfKey;
    using input_type = typename DpfKey::input_type;
    using node_type = typename DpfKey::interior_node;
    using output_type = typename DpfKey::concrete_output_type<I>;

HEDLEY_PRAGMA(GCC diagnostic push)
HEDLEY_PRAGMA(GCC diagnostic ignored "-Wignored-attributes")
    using allocator = aligned_allocator<typename DpfKey::interior_node>;
HEDLEY_PRAGMA(GCC diagnostic pop)
    using unique_ptr = typename allocator::unique_ptr;
    allocator alloc = allocator{};

    if (HEDLEY_UNLIKELY(!std::is_sorted(begin, end)))
    {
        throw std::runtime_error("list must be sorted");
    }
    if (begin == end)
    {
        return subsequence_iterable<DpfKey, decltype(std::begin(outbuf)), ForwardIterator>(
            std::begin(outbuf), begin, end);
    }

    auto mask = dpf_type::msb_mask;
    std::size_t nodes_in_sequence = std::distance(begin, end);
    unique_ptr memo{alloc.allocate_unique_ptr(nodes_in_sequence*2)};

    bool curhalf = (dpf_type::depth ^ 1) & 1;
    memo[!curhalf*nodes_in_sequence + 0] = dpf.root();

    std::vector<ForwardIterator> splits;
    splits.reserve(nodes_in_sequence + 1);
    splits.push_back(begin);
    splits.push_back(end);

    std::size_t level_index = 1;
    auto func = [&](const bool flip = false)
    {
        std::size_t i = 0, j = 0;
        const node_type cw[2] = {
            dpf.correction_word(level_index-1, 0),
            dpf.correction_word(level_index-1, 1)
        };
        const bool is_last = dpf_type::tree::is_last_level(level_index - 1,
            dpf_type::depth);
        // `lower` and `upper` are always adjacent elements of `splits` with `lower` < `upper`
        // [lower, upper) = "block"
        for (std::size_t s = 0; s + 1 < splits.size(); ++s)
        {
            auto lower = splits[s];
            auto upper = splits[s + 1];
            // `upper_bound()` returns iterator to first element where the relevant bit (based on `mask`) is set
            auto it = std::upper_bound(lower, upper, mask,
                [&flip](auto a, auto b){ return static_cast<bool>(a&b) ^ flip; });
            if (it == lower)       // right only since first element in "block" requires right traversal
            {
                memo[curhalf*nodes_in_sequence + i++] = dpf_type::traverse_interior(memo[!curhalf*nodes_in_sequence + j++], cw[1], 1, is_last);
            }
            else if (it == upper)  // left only since no element in "block" requires right traversal
            {
                memo[curhalf*nodes_in_sequence + i++] = dpf_type::traverse_interior(memo[!curhalf*nodes_in_sequence + j++], cw[0], 0, is_last);
            }
            else                    // both ways since some (non-lower) element within "block" requires right traversal
            {
                auto cur_node = memo[!curhalf*nodes_in_sequence + j++];
                auto kids = dpf_type::traverse_interior01(cur_node, cw[0], cw[1], is_last);
                memo[curhalf*nodes_in_sequence + i++] = kids[0];
                memo[curhalf*nodes_in_sequence + i++] = kids[1];
                splits.insert(splits.begin() + static_cast<std::ptrdiff_t>(s + 1), it);
                ++s;  // skip the newly inserted right-block start
            }
        }
    };
    if (dpf_type::depth >= level_index)
    {
        func(utils::uses_signed_msb_v<input_type>);
        ++level_index;
        mask >>= 1;
        curhalf =! curhalf;
    }
    for (; level_index <= dpf_type::depth; ++level_index, mask>>=1, curhalf=!curhalf)
    {
        func();
    }

HEDLEY_PRAGMA(GCC diagnostic push)
HEDLEY_PRAGMA(GCC diagnostic ignored "-Wignored-attributes")
    auto cw = dpf.template leaf<I>();
HEDLEY_PRAGMA(GCC diagnostic pop)
    auto buf = memo.get();

    constexpr auto clz = utils::countl_zero_symmetric_difference<input_type>{};
    auto curr = begin, prev = curr;
    for (std::size_t i = 0, j = 0; i < nodes_in_sequence; ++i)
    {
        j += (clz(*prev, *curr)) < dpf_type::depth;
        auto leaf = dpf_type::template traverse_exterior<I>(buf[j],
            get_if_lo_bit(cw, buf[j]));
        if constexpr (utils::is_packed_subbyte_v<output_type>)
        {
            store_leaf_bytes(outbuf, i, leaf);
        }
        else
        {
            utils::raw_memcpy(&outbuf[i*dpf_type::outputs_per_leaf], &leaf, sizeof(output_type)*dpf_type::outputs_per_leaf);
        }
        prev = curr++;
    }

   return subsequence_iterable<DpfKey, decltype(std::begin(outbuf)), ForwardIterator>(std::begin(outbuf), begin, end);
}

/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          typename DpfKey,
          typename ForwardIterator>
auto eval_sequence_breadth_first(const DpfKey & dpf, ForwardIterator begin, ForwardIterator end)
{
    auto outbuf = make_output_buffer_for_subsequence<I>(dpf, begin, end);

    // moving `outbuf` is allowed as `outbuf` is a `std::vectors`
    //   the underlying data remains on the heap
    //   and thus the data the iterable refers to is still valid
    auto iterable = eval_sequence_breadth_first<I>(dpf, begin, end, outbuf);
    return std::make_pair(std::move(outbuf), std::move(iterable));
}

namespace internal
{

template <typename DpfKey,
          typename SequenceMemoizer>
inline auto eval_sequence_interior(const DpfKey & dpf, const sequence_recipe & recipe,
    SequenceMemoizer && memoizer, std::size_t to_level = DpfKey::depth)
{
    using dpf_type = DpfKey;
    using node_type = typename DpfKey::interior_node;

    // level_index represents the current level being built
    // level_index = 0 => root
    // level_index = depth => last layer of interior nodes
    if (recipe.num_leaf_nodes() == 0)
        return;

    std::size_t level_index = memoizer.assign_dpf(dpf, recipe);
    std::size_t recipe_index = recipe.level_endpoints()[level_index-1];
    std::size_t nodes_at_level = memoizer.get_nodes_at_level(level_index-1);

    for (; level_index <= to_level; level_index = memoizer.advance_level(), nodes_at_level = memoizer.get_nodes_at_level(level_index-1))
    {
        const node_type cw[2] = {
            dpf.correction_word(level_index-1, 0),
            dpf.correction_word(level_index-1, 1)
        };
        const bool is_last = dpf_type::tree::is_last_level(level_index - 1,
            dpf_type::depth);

        auto prevbuf = memoizer[level_index-1];
        auto currbuf = memoizer[level_index];

        DPF_UNROLL_LOOP
        for (std::size_t input_index = 0, output_index = 0; input_index < nodes_at_level; ++input_index, ++recipe_index)
        {
            if (memoizer.traverse_first(recipe_index) == true)
            {
                bool dir = memoizer.get_direction(0);
                currbuf[output_index++] = dpf_type::traverse_interior(prevbuf[input_index], cw[dir], dir, is_last);
            }
            if (memoizer.traverse_second(recipe_index) == true)
            {
                bool dir = memoizer.get_direction(1);
                currbuf[output_index++] = dpf_type::traverse_interior(prevbuf[input_index], cw[dir], dir, is_last);
            }
        }
    }
}

template <std::size_t I,
          typename DpfKey,
          typename OutputBuffer,
          typename SequenceMemoizer>
inline auto eval_sequence_exterior_entire_node(const DpfKey & dpf, const sequence_recipe & recipe,
    OutputBuffer && outbuf, SequenceMemoizer && memoizer)
{
    assert_not_wildcard_output<I>(dpf);

    using dpf_type = DpfKey;
    using output_type = typename DpfKey::concrete_output_type<I>;

    auto nodes_in_interval = recipe.num_leaf_nodes();

HEDLEY_PRAGMA(GCC diagnostic push)
HEDLEY_PRAGMA(GCC diagnostic ignored "-Wignored-attributes")
    auto buf = memoizer[dpf.depth];
HEDLEY_PRAGMA(GCC diagnostic pop)
    DPF_UNROLL_LOOP
    for (std::size_t j = 0; j < nodes_in_interval; ++j)
    {
        auto leaf = dpf.template traverse_exterior<I>(buf[j]);
        if constexpr (utils::is_packed_subbyte_v<output_type>)
        {
            store_leaf_bytes(outbuf, j, leaf);
        }
        else
        {
            utils::raw_memcpy(&outbuf[j*dpf_type::outputs_per_leaf], &leaf, sizeof(output_type)*dpf_type::outputs_per_leaf);
        }
    }
}

template <std::size_t I,
          typename DpfKey,
          typename OutputBuffer,
          typename SequenceMemoizer>
inline auto eval_sequence_exterior_output_only(const DpfKey & dpf, const sequence_recipe & recipe,
    OutputBuffer && outbuf, SequenceMemoizer && memoizer)
{
    assert_not_wildcard_output<I>(dpf);

    using dpf_type = DpfKey;
    using output_type = typename DpfKey::concrete_output_type<I>;

HEDLEY_PRAGMA(GCC diagnostic push)
HEDLEY_PRAGMA(GCC diagnostic ignored "-Wignored-attributes")
    auto cw = dpf.template leaf<I>();
HEDLEY_PRAGMA(GCC diagnostic pop)
    using node_type = typename DpfKey::exterior_node;
    using leaf_node_type = std::tuple_element_t<I, typename DpfKey::leaf_tuple>;
    auto buf = memoizer[dpf.depth];

    leaf_node_type node;
    // DPF_UNROLL_LOOP
    for (std::size_t i = 0, j = -1, prev = -1, curr;
        i < recipe.output_indices().size();
        prev = curr, ++i)
    {
        curr = recipe.output_indices()[i]/dpf_type::outputs_per_leaf;
        if (prev != curr)
        {
            ++j;
            node = dpf_type::template traverse_exterior<I>(buf[j], get_if_lo_bit(cw, buf[j]));
        }
        auto v = extract_leaf<node_type, output_type>(node,
            recipe.output_indices()[i] % dpf_type::outputs_per_leaf);
        assign_share_slot(outbuf[i], v);
    }
}

template <std::size_t ...Is,
          typename DpfKey,
          typename OutputBuffers,
          typename SequenceMemoizer,
          typename ReturnType = return_entire_node_tag_,
          std::enable_if_t<std::is_base_of_v<return_type_tag_, ReturnType>, bool> = true,
          std::size_t ...IIs>
auto eval_sequence(const DpfKey & dpf, const sequence_recipe & recipe,
    OutputBuffers && outbufs, SequenceMemoizer && memoizer, ReturnType return_type, std::index_sequence<IIs...>)
{
    (void)return_type;
    internal::eval_sequence_interior(dpf, recipe, memoizer);

    static_assert(std::is_same_v<ReturnType, return_entire_node_tag_> ||
                    std::is_same_v<ReturnType, return_output_only_tag_>);
    if constexpr (std::is_same_v<ReturnType, return_entire_node_tag_>)
    {
        (internal::eval_sequence_exterior_entire_node<Is>(dpf, recipe, utils::get<IIs>(outbufs), memoizer), ...);
        return utils::make_tuple(
            recipe_subsequence_iterable(std::begin(utils::get<IIs>(outbufs)), recipe.output_indices())...);
    }
    else
    {
        (internal::eval_sequence_exterior_output_only<Is>(dpf, recipe, utils::get<IIs>(outbufs), memoizer), ...);
        const auto nout = recipe.output_indices().size();
        if (nout == 0)
        {
            return utils::make_tuple(subinterval_iterable(std::begin(utils::get<IIs>(outbufs)), utils::size(utils::get<IIs>(outbufs)), 0, 0, 0, 0, false)...);
        }
        return utils::make_tuple(subinterval_iterable(std::begin(utils::get<IIs>(outbufs)), utils::size(utils::get<IIs>(outbufs)), 0, nout-1, 0, 0)...);
    }
}

}  // namespace internal

/// @brief Evaluate `recipe` into a named buffer, reusing `memoizer`.
/// @tparam I output index
/// @tparam Is is
/// @tparam DpfKey DPF key type
/// @tparam OutputBuffers tuple of output buffers
/// @tparam SequenceMemoizer sequence memoizer
/// @tparam ReturnType return type
/// @tparam SequenceMemoizer sequence memoizer
/// @tparam ReturnType return type
/// @param recipe The same object `memoizer` was constructed from.
/// @param outbufs Named buffer. The returned iterable refers into it.
/// @param dpf the DPF key
/// @param memoizer the memoizer built for this key
/// @param return_type `return_entire_node_tag_{}` or `return_output_only_tag_{}`
/// @return the evaluation result
/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          std::size_t ...Is,
          typename DpfKey,
          typename OutputBuffers,
          typename SequenceMemoizer,
          typename ReturnType = return_entire_node_tag_,
          std::enable_if_t<!std::is_base_of_v<return_type_tag_, SequenceMemoizer>, bool> = true,
          std::enable_if_t<std::is_base_of_v<return_type_tag_, ReturnType>, bool> = true>
HEDLEY_ALWAYS_INLINE
auto eval_sequence(const DpfKey & dpf, const sequence_recipe & recipe,
    OutputBuffers & outbufs, SequenceMemoizer && memoizer,  // NOLINT(runtime/references)
    ReturnType return_type = ReturnType{})
{
    assert_not_wildcard_output<I, Is...>(dpf);
    assert_not_wildcard_input(dpf);

    return internal::eval_sequence<I, Is...>(dpf, recipe, outbufs, memoizer, return_type, std::make_index_sequence<1+sizeof...(Is)>());
}

/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          std::size_t ...Is,
          typename DpfKey,
          typename OutputBuffers,
          typename ReturnType = return_entire_node_tag_,
          std::enable_if_t<!std::is_base_of_v<sequence_memoizer_tag_,
              std::decay_t<OutputBuffers>>, bool> = true,
          std::enable_if_t<std::is_base_of_v<return_type_tag_, ReturnType>, bool> = true>
HEDLEY_ALWAYS_INLINE
auto eval_sequence(const DpfKey & dpf, const sequence_recipe & recipe,
    OutputBuffers & outbufs, ReturnType return_type = ReturnType{})  // NOLINT(runtime/references)
{
    return eval_sequence<I, Is...>(dpf, recipe, outbufs,
        dpf::make_double_space_sequence_memoizer<DpfKey>(recipe), return_type);
}

/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          std::size_t ...Is,
          typename DpfKey,
          typename SequenceMemoizer,
          typename ReturnType = return_entire_node_tag_,
          std::enable_if_t<std::is_base_of_v<sequence_memoizer_tag_,
              std::decay_t<SequenceMemoizer>>, bool> = true,
          std::enable_if_t<std::is_base_of_v<return_type_tag_, ReturnType>, bool> = true>
HEDLEY_ALWAYS_INLINE
auto eval_sequence(const DpfKey & dpf, const sequence_recipe & recipe,
    SequenceMemoizer && memoizer, ReturnType return_type = ReturnType{})
{
    auto outbufs = utils::make_tuple(
        make_output_buffer_for_recipe_subsequence<I>(dpf, recipe, return_type),
        make_output_buffer_for_recipe_subsequence<Is>(dpf, recipe, return_type)...);

    // moving `outbufs` is allowed as the `outbufs` are `std::vectors`
    //   the underlying data remains on the heap
    //   and thus the data the iterable refers to is still valid
    auto iterable = eval_sequence<I, Is...>(dpf, recipe, outbufs, memoizer, return_type);
    return std::make_pair(std::move(outbufs), std::move(iterable));
}

/// \complexity O(n k) interior traversals in the worst case and O(k) node workspace. k is the number of listed points and n is `depth`. The breadth-first buffer is 2k nodes, so each level traverses at most one node per point. Shared prefixes do fewer traversals. A recipe memoizer instead stores O(recipe leaf nodes) (see that memoizer).
template <std::size_t I = 0,
          std::size_t ...Is,
          typename DpfKey,
          typename ReturnType = return_entire_node_tag_,
          std::enable_if_t<std::is_base_of_v<return_type_tag_, ReturnType>, bool> = true>
HEDLEY_ALWAYS_INLINE
auto eval_sequence(const DpfKey & dpf, const sequence_recipe & recipe,
     ReturnType return_type = ReturnType{})
{
    return eval_sequence<I, Is...>(dpf, recipe,
       dpf::make_double_space_sequence_memoizer<DpfKey>(recipe), return_type);
}

/// @brief Fold a sorted sequence into `pr` via interval-run covers.
/// @details Maximal contiguous runs use once-per-BFS-node `prove_fold_interval`.
///          Isolated points are length-1 runs. Both parties must see the same
///          sorted list. Prefer this over path-memo when the list is dense.
template <typename KeyT, typename ForwardIterator>
void prove_sequence(const KeyT & key, ForwardIterator begin, ForwardIterator end,
    prove_ref pr)
{
    static_assert(KeyT::is_verifiable,
        "prove_sequence: key must carry dpf::verifiable");
    if (HEDLEY_UNLIKELY(begin != end && !std::is_sorted(begin, end)))
        throw std::runtime_error("list must be sorted");
    detail::vdpf::init_proof(pr.token, key);
    using input_type = typename KeyT::input_type;
    for (auto it = begin; it != end; )
    {
        const auto run_from = static_cast<input_type>(*it);
        auto run_to = run_from;
        ++it;
        while (it != end)
        {
            const auto next = static_cast<input_type>(*it);
            if (next != static_cast<input_type>(run_to + input_type{1}))
                break;
            run_to = next;
            ++it;
        }
        prove_fold_interval(key, run_from, run_to, pr.token);
    }
    detail::vdpf::fold_output_binding(pr.token, key);
}

}  // namespace dpf

#endif  // LIBDPF_INCLUDE_DPF_EVAL_SEQUENCE_HPP__