boost::histogram::axis::variant, allow users to choose between sorted_array+std::lower_bound and eytzinger_layout+eytzinger_binary_search

Inizia individuando boost::histogram::axis::variant e il test a cui si fa riferimento nell’issue, quindi esamina come viene attualmente utilizzato std::upper_bound. Confronta i due layout di ricerca proposti e determina come gli utenti sceglierebbero tra loro; il lavoro è completato quando la scelta è esposta senza modificare il comportamento esistente e i test o benchmark pertinenti coprono entrambe le opzioni.

According to the test,

• When input follows normal distribution, sorted_array+std::lower_bound is faster;
• When input follows uniform distribution, eytzinger_layout+eytzinger_binary_search is faster;
  Eytzinger Binary Search - Algorithmica
  Static B-Trees: up to 15x faster than std::lower_bound : cpp

(boost::histogram::axis::variant is currently using std::upper_bound.)

comparison (x86-64 gcc (trunk), -std=c++23 -O3 -Wall -pedantic -pthread):

#pragma GCC optimize("O3")
#include <bits/stdc++.h>

constexpr size_t n = (1<<20);
constexpr size_t block_size = 64 / sizeof(int); // = 64 / 4 = cache_line_size / sizeof(int)
alignas(64) int a[n], b[n+1];

constexpr size_t query_count=(1<<20);
int targets[query_count];
int results_eytzinger[query_count];
int results_lower_bound[query_count];

int eytzinger(int i = 0, size_t k = 1) {
    if (k <= n) {
        i = eytzinger(i, 2 * k);
        b[k] = a[i++];
        i = eytzinger(i, 2 * k + 1);
    }
    return i;
}

int search(int x) {
    size_t k = 1;
    while (k <= n) {
        __builtin_prefetch(b + k * block_size);
        k = 2 * k + (b[k] < x);
    }
    k >>= __builtin_ffs(~k);
    return k;
}

int main()
{
    std::random_device rd;  //Will be used to obtain a seed for the random number engine
    std::mt19937 gen(rd()); //Standard mersenne_twister_engine seeded with rd()
    std::uniform_int_distribution<> distrib(std::numeric_limits<int>::min(), std::numeric_limits<int>::max());
    for (size_t i = 0; i != n; ++i)
    {
        a[i]=std::round(distrib(gen));
    }
    std::sort(std::begin(a),std::end(a));

    eytzinger();

    for (size_t i = 0; i != query_count; ++i)
    {
        targets[i]=std::round(distrib(gen));
    }

    {
        auto start = std::chrono::steady_clock::now();
        for (size_t i = 0; i != query_count; ++i)
        {
            results_eytzinger[i]=search(targets[i]);
        }
        auto end = std::chrono::steady_clock::now();
        std::chrono::duration<double> elapsed_seconds = end-start;
        std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
    }
    {
        auto start = std::chrono::steady_clock::now();
        for (size_t i = 0; i != query_count; ++i)
        {
            results_lower_bound[i]=std::lower_bound(std::begin(a),std::end(a),targets[i])-std::begin(a);
        }
        auto end = std::chrono::steady_clock::now();
        std::chrono::duration<double> elapsed_seconds = end-start;
        std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
    }
    for (size_t i = 0; i != query_count; ++i)
    {
        assert(a[results_lower_bound[i]]==b[results_eytzinger[i]]);   
    }
}