add eytzinger to benchmarks

benchmark/CMakeLists.txt

@@ -8,6 +8,7 @@ include(BoostFetch)
 set(CMAKE_BUILD_TYPE Release) # ok, only set in local scope
 option(BENCHMARK_ENABLE_TESTING "" OFF)
 boost_fetch(google/benchmark TAG main)
+boost_fetch(boostorg/align TAG develop EXCLUDE_FROM_ALL)
 boost_fetch(boostorg/math TAG develop EXCLUDE_FROM_ALL)
 
 function(add_benchmark NAME)
@@ -27,10 +28,10 @@ function(add_benchmark NAME)
 
   add_executable(${NAME} ${SOURCE})
   target_include_directories(${NAME} PRIVATE ${__INCLUDE_DIRECTORIES})
-  target_link_libraries(${NAME} PRIVATE Boost::histogram Boost::math benchmark_main ${__LINK_LIBRARIES})
-  target_compile_options(${NAME} PRIVATE -DNDEBUG -O3 -funsafe-math-optimizations ${__COMPILE_OPTIONS})
-  if (NOT DARWIN)
-    target_compile_options(${NAME} PRIVATE -march=native)
+  target_link_libraries(${NAME} PRIVATE Boost::histogram Boost::math Boost::align benchmark_main ${__LINK_LIBRARIES})
+  target_compile_options(${NAME} PRIVATE -DNDEBUG -O3 -funsafe-math-optimizations ${__COMPILE_OPTIONS} -DHAVE_BOOST_ALIGN)
+  if (NOT CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
+     target_compile_options(${NAME} PRIVATE -march=native ${__COMPILE_OPTIONS})
   endif()
 endfunction()
 
@@ -39,6 +40,7 @@ add_benchmark(axis_index)
 add_benchmark(histogram_filling)
 add_benchmark(histogram_iteration)
 add_benchmark(detail_normal)
+add_benchmark(detail_eytzinger_search)
 
 find_package(Threads)
 if (Threads_FOUND)

benchmark/detail_eytzinger_search.cpp

@@ -0,0 +1,116 @@
+// Copyright 2022 Hans Dembinski
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <benchmark/benchmark.h>
+#include <boost/histogram/axis/variable.hpp>
+#include <numeric>
+#include "../test/throw_exception.hpp"
+#include "generator.hpp"
+#include <algorithm>
+#include <vector>
+#include <iostream>
+#include <boost/core/bit.hpp>
+#ifdef HAVE_BOOST_ALIGN
+  #include <boost/align/aligned_allocator.hpp>
+#endif
+#include <boost/core/lightweight_test.hpp>
+
+#include <cassert>
+struct assert_check {
+  assert_check() {
+    assert(false); // don't run with asserts enabled
+  }
+} _;
+
+template<typename T=double>
+struct eytzinger_search {
+
+  int ffs(size_t v) const noexcept
+  {
+      if(v==0) return 0;
+  #if HAVE_BOOST_MULTIPRECISION
+      return boost::multiprecision::lsb(v)+1;
+  #else
+      // we prefer boost::core since it is a dependency already
+      return boost::core::countr_zero(v)+1;
+  #endif
+  }
+
+  eytzinger_search(const std::vector<T>& a) :
+      b_(a.size() + 1), idx_(a.size() + 1)
+  {
+      init(a);
+      idx_[0] = a.size() - 1;
+  }
+
+  int index(T const& x) const {
+      size_t k = 1;
+      while (k < b_.size())
+          k = 2 * k + (b_[k] < x);
+      k >>= ffs(~k);
+      return idx_[k];
+  }
+
+  size_t init(const std::vector<T>& a, size_t i = 0, size_t k = 1) {
+      if (k <= a.size()) {
+          i = init(a, i, 2 * k);
+          idx_[k] = i - 1;
+          b_[k] = a[i++];
+          i = init(a, i, 2 * k + 1);
+      }
+      return i;
+  }
+
+#ifdef HAVE_BOOST_ALIGN
+  std::vector<T, boost::alignment::aligned_allocator<T, 64>> b_;
+  std::vector<int, boost::alignment::aligned_allocator<int, 64>> idx_;
+#else
+  std::vector<T> b_;
+  std::vector<int> idx_;
+#endif
+};
+
+using namespace boost::histogram;
+
+template <class Distribution>
+static void variable(benchmark::State& state) {
+  std::vector<double> v;
+  for (double x = 0; x <= state.range(0); ++x) { v.push_back(x / state.range(0)); }
+  auto a = axis::variable<>(v);
+  generator<Distribution> gen;  
+  for (auto _ : state) benchmark::DoNotOptimize(a.index(gen()));
+}
+
+template <class Distribution>
+static void eytzinger(benchmark::State& state) {
+  std::vector<double> v;
+  for (double x = 0; x <= state.range(0); ++x) { v.push_back(x / state.range(0)); }
+  auto a = eytzinger_search<double>(v);
+  generator<Distribution> gen;    
+
+  auto aref = axis::variable<>(v);
+  for (int i = 0; i < 10000; ++i) {
+    const double x = gen();
+    if (a.index(x) != aref.index(x)) {
+      std::cout << "x = " << x << " " 
+                << a.index(x) << " " 
+                << aref.index(x) << std::endl;
+      std::abort();
+    }
+  }
+
+  for (auto _ : state) benchmark::DoNotOptimize(a.index(gen()));
+}
+
+BENCHMARK_TEMPLATE(variable, uniform)->RangeMultiplier(10)->Range(10, 10000);
+BENCHMARK_TEMPLATE(variable, normal)->RangeMultiplier(10)->Range(10, 10000);
+BENCHMARK_TEMPLATE(variable, chi2)->RangeMultiplier(10)->Range(10, 10000);
+BENCHMARK_TEMPLATE(variable, expon)->RangeMultiplier(10)->Range(10, 10000);
+
+BENCHMARK_TEMPLATE(eytzinger, uniform)->RangeMultiplier(10)->Range(10, 10000);
+BENCHMARK_TEMPLATE(eytzinger, normal)->RangeMultiplier(10)->Range(10, 10000);
+BENCHMARK_TEMPLATE(eytzinger, chi2)->RangeMultiplier(10)->Range(10, 10000);
+BENCHMARK_TEMPLATE(eytzinger, expon)->RangeMultiplier(10)->Range(10, 10000);

benchmark/generator.hpp

@@ -10,6 +10,8 @@
 using uniform = std::uniform_real_distribution<>;
 using uniform_int = std::uniform_int_distribution<>;
 using normal = std::normal_distribution<>;
+using chi2 = std::chi_squared_distribution<>;
+using expon = std::exponential_distribution<>;
 
 template <class Distribution, class... Ts>
 Distribution init(Ts...);
@@ -24,6 +26,16 @@ normal init<normal>() {
   return normal{0.5, 0.3};
 }
 
+template <>
+chi2 init<chi2>() {
+  return chi2{5.0};
+}
+
+template <>
+expon init<expon>() {
+  return expon{1};
+}
+
 template <>
 uniform_int init<uniform_int, int>(int n) {
   return uniform_int{0, n};