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,122 @@
+// 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 <algorithm>
+#include <boost/core/bit.hpp>
+#include <boost/histogram/axis/variable.hpp>
+#include <boost/histogram/detail/prefetch.hpp>
+#include <iostream>
+#include <numeric>
+#include <vector>
+#include "../test/throw_exception.hpp"
+#include "generator.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
+  }
+} _;
+
+using namespace boost::histogram;
+
+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()) {
+      constexpr int block_size = detail::cacheline_length / sizeof(T);
+      detail::prefetch(b_.data() + k * block_size);
+      k = 2 * k + !(x < b_[k]); // double negation to handle nan correctly
+    }
+    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, detail::cache_alignment>> b_;
+  std::vector<int, boost::alignment::aligned_allocator<int, detail::cache_alignment>>
+      idx_;
+#else
+  std::vector<T> b_;
+  std::vector<int> idx_;
+#endif
+};
+
+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);
+  std::vector<double> inputs(1000);
+  std::generate(inputs.begin(), inputs.end(), gen);
+  inputs[0] = std::numeric_limits<double>::quiet_NaN();
+  inputs[1] = -std::numeric_limits<double>::infinity();
+  inputs[2] = std::numeric_limits<double>::infinity();
+  inputs[3] = v[0];
+  inputs[4] = v.back();
+  for (auto&& x : inputs) {
+    if (a.index(x) != aref.index(x)) {
+      std::cout << "x = " << x << ": new = " << a.index(x) << " ref = " << aref.index(x)
+                << std::endl;
+      std::abort();
+    }
+  }
+
+  for (auto _ : state) benchmark::DoNotOptimize(a.index(gen()));
+}
+
+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()));
+}
+
+BENCHMARK_TEMPLATE(eytzinger, uniform)->RangeMultiplier(10)->Range(10, 100000);
+BENCHMARK_TEMPLATE(eytzinger, normal)->RangeMultiplier(10)->Range(10, 100000);
+BENCHMARK_TEMPLATE(eytzinger, chi2)->RangeMultiplier(10)->Range(10, 100000);
+BENCHMARK_TEMPLATE(eytzinger, expon)->RangeMultiplier(10)->Range(10, 100000);
+
+BENCHMARK_TEMPLATE(variable, uniform)->RangeMultiplier(10)->Range(10, 100000);
+BENCHMARK_TEMPLATE(variable, normal)->RangeMultiplier(10)->Range(10, 100000);
+BENCHMARK_TEMPLATE(variable, chi2)->RangeMultiplier(10)->Range(10, 100000);
+BENCHMARK_TEMPLATE(variable, expon)->RangeMultiplier(10)->Range(10, 100000);

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};

include/boost/histogram/detail/prefetch.hpp

@@ -0,0 +1,54 @@
+//          Copyright Oliver Kowalke, Hans Dembinski 2014-2022.
+// 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)
+
+#ifndef BOOST_HISTOGRAM_DETAIL_PREFETCH_H
+#define BOOST_HISTOGRAM_DETAIL_PREFETCH_H
+
+#include <boost/config.hpp>
+#include <cstddef>
+#include <cstdint>
+
+#if BOOST_GCC || BOOST_CLANG
+#define BOOST_HISTOGRAM_DETAIL_HAVE_BUILTIN_PREFETCH
+#elif BOOST_INTEL
+#define BOOST_HISTOGRAM_DETAIL_HAVE_MM_PREFETCH
+#include <immintrin.h>
+#elif BOOST_MSVC && !defined(_M_ARM) && !defined(_M_ARM64)
+#define BOOST_HISTOGRAM_DETAIL_HAVE_MM_PREFETCH
+#include <immintrin.h>
+#endif
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+// modern architectures have cachelines with 64byte length
+// ARM Cortex-A15 32/64byte, Cortex-A9 16/32/64bytes
+// MIPS 74K: 32byte, 4KEc: 16byte
+// ist should be safe to use 64byte for all
+static constexpr std::size_t cache_alignment{64};
+static constexpr std::size_t cacheline_length{64};
+// lookahead size for prefetching
+static constexpr std::size_t prefetch_stride{4 * cacheline_length};
+
+BOOST_FORCEINLINE
+void prefetch(const void* addr) {
+#if defined BOOST_HISTOGRAM_DETAIL_HAVE_MM_PREFETCH
+  _mm_prefetch((const char*)addr, _MM_HINT_T0);
+#elif defined BOOST_HISTOGRAM_DETAIL_HAVE_BUILTIN_PREFETCH
+  // L1 cache : hint == 1
+  __builtin_prefetch((void*)addr, 1, 1);
+#endif
+  // no prefetch available
+}
+
+#undef BOOST_HISTOGRAM_DETAIL_HAVE_MM_PREFETCH
+#undef BOOST_HISTOGRAM_DETAIL_HAVE_BUILTIN_PREFETCH
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif // BOOST_HISTOGRAM_DETAIL_PREFETCH_H

test/odr_test.cpp

@@ -17,6 +17,7 @@
 
 #include <boost/histogram/detail/ignore_deprecation_warning_begin.hpp>
 #include <boost/histogram/detail/ignore_deprecation_warning_end.hpp>
+#include <boost/histogram/detail/prefetch.hpp>
 #include <boost/histogram/utility/clopper_pearson_interval.hpp>
 #include <boost/histogram/utility/jeffreys_interval.hpp>
 #include <boost/histogram/utility/wald_interval.hpp>