From d5dd9ec39bb0e3e266ab8bef5c6f1e2c1474bddd Mon Sep 17 00:00:00 2001
From: Hans Dembinski <hans.dembinski@gmail.com>
Date: Thu, 22 Dec 2022 13:35:31 +0100
Subject: [PATCH 1/5] add eytzinger to benchmarks

---
 benchmark/CMakeLists.txt              |  10 ++-
 benchmark/detail_eytzinger_search.cpp | 116 ++++++++++++++++++++++++++
 benchmark/generator.hpp               |  12 +++
 3 files changed, 134 insertions(+), 4 deletions(-)
 create mode 100644 benchmark/detail_eytzinger_search.cpp

diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt
index 656e6f08d..ed4e7ecba 100644
--- a/benchmark/CMakeLists.txt
+++ b/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)
diff --git a/benchmark/detail_eytzinger_search.cpp b/benchmark/detail_eytzinger_search.cpp
new file mode 100644
index 000000000..5a37614c7
--- /dev/null
+++ b/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);
diff --git a/benchmark/generator.hpp b/benchmark/generator.hpp
index f613ced35..3cacc3224 100644
--- a/benchmark/generator.hpp
+++ b/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};

From 8ab311f000c2c0e363ebecf5db18495707b25853 Mon Sep 17 00:00:00 2001
From: Hans Dembinski <hans.dembinski@gmail.com>
Date: Fri, 23 Dec 2022 11:30:34 +0100
Subject: [PATCH 2/5] support prefetch

---
 benchmark/detail_eytzinger_search.cpp       | 120 +++++++++++---------
 include/boost/histogram/detail/prefetch.hpp |  54 +++++++++
 2 files changed, 120 insertions(+), 54 deletions(-)
 create mode 100644 include/boost/histogram/detail/prefetch.hpp

diff --git a/benchmark/detail_eytzinger_search.cpp b/benchmark/detail_eytzinger_search.cpp
index 5a37614c7..d4f8d1823 100644
--- a/benchmark/detail_eytzinger_search.cpp
+++ b/benchmark/detail_eytzinger_search.cpp
@@ -5,16 +5,17 @@
 // 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"
-#include <algorithm>
-#include <vector>
-#include <iostream>
-#include <boost/core/bit.hpp>
 #ifdef HAVE_BOOST_ALIGN
-  #include <boost/align/aligned_allocator.hpp>
+#include <boost/align/aligned_allocator.hpp>
 #endif
 #include <boost/core/lightweight_test.hpp>
 
@@ -25,63 +26,66 @@ struct assert_check {
   }
 } _;
 
-template<typename T=double>
-struct eytzinger_search {
+using namespace boost::histogram;
 
-  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
+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;
+  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 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;
+    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_;
+  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
 };
 
-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()));
+  generator<Distribution> gen;
+  std::vector<double> x(10000);
+  std::generate(x.begin(), x.end(), gen);
+  unsigned i = 0;
+  for (auto _ : state) benchmark::DoNotOptimize(a.index(x[i++ % x.size()]));
 }
 
 template <class Distribution>
@@ -89,28 +93,36 @@ 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;    
+  generator<Distribution> gen;
 
   auto aref = axis::variable<>(v);
-  for (int i = 0; i < 10000; ++i) {
-    const double x = gen();
+  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 << " " 
-                << a.index(x) << " " 
-                << aref.index(x) << std::endl;
+      std::cout << "x = " << x << ": new = " << a.index(x) << " ref = " << aref.index(x)
+                << std::endl;
       std::abort();
     }
   }
 
-  for (auto _ : state) benchmark::DoNotOptimize(a.index(gen()));
+  std::vector<double> x(10000);
+  std::generate(x.begin(), x.end(), gen);
+  unsigned i = 0;
+  for (auto _ : state) benchmark::DoNotOptimize(a.index(x[i++ % x.size()]));
 }
 
-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, 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(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_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);
diff --git a/include/boost/histogram/detail/prefetch.hpp b/include/boost/histogram/detail/prefetch.hpp
new file mode 100644
index 000000000..8fde26ab3
--- /dev/null
+++ b/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 <mmintrin.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

From df791a4bd95432fda5ee50a83dc230d6cd981585 Mon Sep 17 00:00:00 2001
From: Hans Dembinski <hans.dembinski@gmail.com>
Date: Fri, 23 Dec 2022 11:33:21 +0100
Subject: [PATCH 3/5] back to generate values live

---
 benchmark/detail_eytzinger_search.cpp | 26 ++++++++++----------------
 1 file changed, 10 insertions(+), 16 deletions(-)

diff --git a/benchmark/detail_eytzinger_search.cpp b/benchmark/detail_eytzinger_search.cpp
index d4f8d1823..ff00bf37c 100644
--- a/benchmark/detail_eytzinger_search.cpp
+++ b/benchmark/detail_eytzinger_search.cpp
@@ -76,18 +76,6 @@ struct eytzinger_search {
 #endif
 };
 
-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;
-  std::vector<double> x(10000);
-  std::generate(x.begin(), x.end(), gen);
-  unsigned i = 0;
-  for (auto _ : state) benchmark::DoNotOptimize(a.index(x[i++ % x.size()]));
-}
-
 template <class Distribution>
 static void eytzinger(benchmark::State& state) {
   std::vector<double> v;
@@ -111,10 +99,16 @@ static void eytzinger(benchmark::State& state) {
     }
   }
 
-  std::vector<double> x(10000);
-  std::generate(x.begin(), x.end(), gen);
-  unsigned i = 0;
-  for (auto _ : state) benchmark::DoNotOptimize(a.index(x[i++ % x.size()]));
+  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);

From fda62a198884d902367fe0317197dcf4c358facc Mon Sep 17 00:00:00 2001
From: Hans Dembinski <hans.dembinski@gmail.com>
Date: Fri, 23 Dec 2022 11:43:28 +0100
Subject: [PATCH 4/5] fix

---
 test/odr_test.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/odr_test.cpp b/test/odr_test.cpp
index 13b108bc0..435ecd186 100644
--- a/test/odr_test.cpp
+++ b/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>

From 2caa72498239a8e6a6e31b487349805c8b2ab9f0 Mon Sep 17 00:00:00 2001
From: Hans Dembinski <hans.dembinski@gmail.com>
Date: Sat, 24 Dec 2022 01:29:41 +0100
Subject: [PATCH 5/5] try to fix windows CI

---
 include/boost/histogram/detail/prefetch.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/boost/histogram/detail/prefetch.hpp b/include/boost/histogram/detail/prefetch.hpp
index 8fde26ab3..470a46cd9 100644
--- a/include/boost/histogram/detail/prefetch.hpp
+++ b/include/boost/histogram/detail/prefetch.hpp
@@ -17,7 +17,7 @@
 #include <immintrin.h>
 #elif BOOST_MSVC && !defined(_M_ARM) && !defined(_M_ARM64)
 #define BOOST_HISTOGRAM_DETAIL_HAVE_MM_PREFETCH
-#include <mmintrin.h>
+#include <immintrin.h>
 #endif
 
 namespace boost {