Fractional Cascading: Searching K Sorted Lists in O(log N + K)

The $O(k \log n) \to O(\log n + k)$ Iterative Search Bottleneck

In computational geometry (2D orthogonal range trees, segment trees of sorted vectors), database multi-column index lookups, and packet classification, we often store a catalog of $k$ sorted arrays $L_0, L_1, \dots, L_{k-1}$, each of size $\approx n$, and need to compute $\text{lower\_bound}(L_i, q)$ for the same query value $q$ in every list $i = 0 \dots k-1$.

Executing $k$ independent binary searches takes $O(k \log n)$ time. Could we store a pointer from each element of $L_0$ to its successor in $L_1$, $L_2$, etc.? If we simply merged all $k$ lists into $L_0$, $L_0$ would blow up to size $k n$, and propagating full lists across $k$ levels would consume $\Theta(k^2 n)$ space! Bernard Chazelle and Leonidas Guibas (1986) solved this space-time trade-off with Fractional Cascading: propagate only a $\frac{1}{2}$-fraction of each augmented list upward!

Production C++20 Implementation: Complete FractionalCascadingCatalog

#include <algorithm>
#include <concepts>
#include <cstddef>
#include <limits>
#include <vector>

template <std::totally_ordered T>
class FractionalCascadingCatalog {
public:
    struct Node {
        T val;
        std::size_t real_idx;         // Index of lower_bound(val) in original list L_i
        std::size_t next_cascade_idx; // Index of lower_bound(val) in next augmented list M_{i+1}
    };

    explicit FractionalCascadingCatalog(const std::vector<std::vector<T>>& lists, T sentinel_max)
        : k_(lists.size()), sentinel_(sentinel_max), orig_(lists), aug_(k_) {
        if (k_ == 0) return;

        // Build bottom-up from level k - 1 down to 0
        for (std::ptrdiff_t lvl = static_cast<std::ptrdiff_t>(k_) - 1; lvl >= 0; --lvl) {
            const auto& L = orig_[lvl];
            std::vector<T> merged;
            if (static_cast<std::size_t>(lvl) == k_ - 1) {
                merged = L;
            } else {
                // Extract every 2nd element (indices 1, 3, 5, ...) of M_{lvl+1}
                std::vector<T> promoted;
                const auto& M_next = aug_[lvl + 1];
                for (std::size_t j = 1; j < M_next.size(); j += 2) {
                    if (M_next[j].val < sentinel_) promoted.push_back(M_next[j].val);
                }
                merged.resize(L.size() + promoted.size());
                std::merge(L.begin(), L.end(), promoted.begin(), promoted.end(), merged.begin());
                merged.erase(std::unique(merged.begin(), merged.end()), merged.end());
            }
            merged.push_back(sentinel_); // +infinity sentinel at end

            // Populate real_idx (into L) and next_cascade_idx (into M_{lvl+1}) via two-pointer sweep
            auto& M = aug_[lvl];
            M.resize(merged.size());
            std::size_t r_ptr = 0, c_ptr = 0;
            for (std::size_t idx = 0; idx < merged.size(); ++idx) {
                const T v = merged[idx];
                while (r_ptr < L.size() && L[r_ptr] < v) ++r_ptr;
                if (static_cast<std::size_t>(lvl) + 1 < k_) {
                    const auto& M_next = aug_[lvl + 1];
                    while (c_ptr + 1 < M_next.size() && M_next[c_ptr].val < v) ++c_ptr;
                }
                M[idx] = Node{v, r_ptr, c_ptr};
            }
        }
    }

    // Returns the lower_bound index in each of the k original lists in O(log N + k) time!
    std::vector<std::size_t> query_all_lower_bounds(const T& q) const {
        if (k_ == 0) return {};
        std::vector<std::size_t> results(k_);

        // 1. Single O(log N) binary search in top augmented list M_0
        const auto& M0 = aug_[0];
        auto it = std::lower_bound(M0.begin(), M0.end(), q,
                                   [](const Node& nd, const T& val) { return nd.val < val; });
        std::size_t cur = static_cast<std::size_t>(it - M0.begin());
        results[0] = M0[cur].real_idx;

        // 2. O(1) cascade step for each subsequent list 1 .. k - 1
        for (std::size_t lvl = 0; lvl + 1 < k_; ++lvl) {
            std::size_t next_pos = aug_[lvl][cur].next_cascade_idx;
            const auto& M_next = aug_[lvl + 1];
            // Because every 2nd element of M_{lvl+1} was promoted to M_lvl,
            // the exact lower_bound in M_{lvl+1} is either next_pos - 1 or next_pos!
            if (next_pos > 0 && !(M_next[next_pos - 1].val < q)) {
                --next_pos;
            }
            cur = next_pos;
            results[lvl + 1] = M_next[cur].real_idx;
        }
        return results;
    }

private:
    std::size_t k_;
    T sentinel_;
    std::vector<std::vector<T>> orig_;
    std::vector<std::vector<Node>> aug_;
};