diff --git a/CMakeLists.txt b/CMakeLists.txt
index 78196d49..960b7252 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -94,6 +94,7 @@ compiled_test(test/axis_variant_test.cpp)
 compiled_test(test/detail_test.cpp)
 compiled_test(test/histogram_dynamic_test.cpp)
 compiled_test(test/histogram_mixed_test.cpp)
+compiled_test(test/histogram_operators_test.cpp)
 compiled_test(test/histogram_test.cpp)
 compiled_test(test/indexed_test.cpp)
 compiled_test(test/internal_accumulators_test.cpp)
diff --git a/examples/guide_histogram_operators.cpp b/examples/guide_histogram_operators.cpp
index 6d2d5d6c..b80d1937 100644
--- a/examples/guide_histogram_operators.cpp
+++ b/examples/guide_histogram_operators.cpp
@@ -24,23 +24,24 @@ int main() {
   auto h3 = h1;
   h3 += h2; // counts are: 1 1
 
-  // adding multiple histograms at once is efficient and does not create
-  // superfluous temporaries since operator+ functions are overloaded to
-  // accept and return rvalue references where possible
+  // adding multiple histograms at once is likely to be optimized by the compiler so that
+  // superfluous temporaries avoided, but no guarantees are given; use this equivalent
+  // code when you want to make sure: h4 = h1; h4 += h2; h4 += h3;
   auto h4 = h1 + h2 + h3; // counts are: 2 2
 
   assert(h4.at(0) == 2 && h4.at(1) == 2);
 
-  // multiply by number
-  h4 *= 2; // counts are: 4 4
+  // multiply by number; h4 *= 2 is not allowed, because the result of a multiplication is
+  // not a histogram anymore, it has the type boost::histogram::grid<...>)
+  auto g4 = h4 * 2; // counts are: 4 4
 
-  // divide by number
-  auto h5 = h4 / 4; // counts are: 1 1
+  // divide by number; g4 /= 4 also works
+  auto g5 = g4 / 4; // counts are: 1 1
 
-  assert(h5.at(0) == 1 && h5.at(1) == 1);
-  assert(h4 != h5 && h4 == 4 * h5);
+  assert(g5.at(0) == 1 && g5.at(1) == 1);
+  assert(g4 != g5 && g4 == 4 * g5);
 
-  // note special effect of multiplication on weighted_sum variance
+  // note special effect of multiplication on weighted_sum
   auto h = bh::make_histogram_with(std::vector<bh::accumulators::weighted_sum<double>>(),
                                    bh::axis::regular<>(2, -1, 1));
   h(-0.5);
@@ -49,14 +50,14 @@ int main() {
   assert(h.at(0).value() == 1 && h.at(1).value() == 0);
 
   auto h_sum = h + h;
-  auto h_mul = 2 * h;
+  auto g_mul = 2 * h;
 
   // equality operator checks variances, so following statement is false
-  assert(h_sum != h_mul);
+  assert(h_sum != g_mul);
 
   // variance is different when histograms are scaled
-  assert(h_sum.at(0).value() == 2 && h_mul.at(0).value() == 2);
-  assert(h_sum.at(0).variance() == 2 && h_mul.at(0).variance() == 4);
+  assert(h_sum.at(0).value() == 2 && g_mul.at(0).value() == 2);
+  assert(h_sum.at(0).variance() == 2 && g_mul.at(0).variance() == 4);
 }
 
 //]
diff --git a/include/boost/histogram.hpp b/include/boost/histogram.hpp
index 8b7f693f..fccfc0b9 100644
--- a/include/boost/histogram.hpp
+++ b/include/boost/histogram.hpp
@@ -16,6 +16,7 @@
 #include <boost/histogram/algorithm/reduce.hpp>
 #include <boost/histogram/algorithm/sum.hpp>
 #include <boost/histogram/axis.hpp>
+#include <boost/histogram/grid.hpp>
 #include <boost/histogram/histogram.hpp>
 #include <boost/histogram/indexed.hpp>
 #include <boost/histogram/literals.hpp>
diff --git a/include/boost/histogram/accumulators/mean.hpp b/include/boost/histogram/accumulators/mean.hpp
index f811086d..540c19ae 100644
--- a/include/boost/histogram/accumulators/mean.hpp
+++ b/include/boost/histogram/accumulators/mean.hpp
@@ -7,8 +7,9 @@
 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_MEAN_HPP
 #define BOOST_HISTOGRAM_ACCUMULATORS_MEAN_HPP
 
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <cstddef>
+#include <type_traits>
 
 namespace boost {
 namespace histogram {
@@ -75,4 +76,12 @@ private:
 } // namespace histogram
 } // namespace boost
 
+namespace std {
+template <class T, class U>
+struct common_type<boost::histogram::accumulators::mean<T>,
+                   boost::histogram::accumulators::mean<U>> {
+  using type = boost::histogram::accumulators::mean<common_type_t<T, U>>;
+};
+} // namespace std
+
 #endif
diff --git a/include/boost/histogram/accumulators/ostream_operators.hpp b/include/boost/histogram/accumulators/ostream_operators.hpp
index 899ebd8c..79982acd 100644
--- a/include/boost/histogram/accumulators/ostream_operators.hpp
+++ b/include/boost/histogram/accumulators/ostream_operators.hpp
@@ -7,7 +7,7 @@
 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_OSTREAM_OPERATORS_HPP
 #define BOOST_HISTOGRAM_ACCUMULATORS_OSTREAM_OPERATORS_HPP
 
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <ostream>
 
 namespace boost {
diff --git a/include/boost/histogram/accumulators/sum.hpp b/include/boost/histogram/accumulators/sum.hpp
index 34ad4bb9..cc4fc71f 100644
--- a/include/boost/histogram/accumulators/sum.hpp
+++ b/include/boost/histogram/accumulators/sum.hpp
@@ -7,8 +7,9 @@
 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
 #define BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
 
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <cmath>
+#include <type_traits>
 
 namespace boost {
 namespace histogram {
@@ -84,4 +85,12 @@ private:
 } // namespace histogram
 } // namespace boost
 
+namespace std {
+template <class T, class U>
+struct common_type<boost::histogram::accumulators::sum<T>,
+                   boost::histogram::accumulators::sum<U>> {
+  using type = boost::histogram::accumulators::sum<common_type_t<T, U>>;
+};
+} // namespace std
+
 #endif
diff --git a/include/boost/histogram/accumulators/weighted_mean.hpp b/include/boost/histogram/accumulators/weighted_mean.hpp
index 9bd8be3c..574058c8 100644
--- a/include/boost/histogram/accumulators/weighted_mean.hpp
+++ b/include/boost/histogram/accumulators/weighted_mean.hpp
@@ -7,7 +7,8 @@
 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_MEAN_HPP
 #define BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_MEAN_HPP
 
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
+#include <type_traits>
 
 namespace boost {
 namespace histogram {
@@ -79,4 +80,12 @@ private:
 } // namespace histogram
 } // namespace boost
 
+namespace std {
+template <class T, class U>
+struct common_type<boost::histogram::accumulators::weighted_mean<T>,
+                   boost::histogram::accumulators::weighted_mean<U>> {
+  using type = boost::histogram::accumulators::weighted_mean<common_type_t<T, U>>;
+};
+} // namespace std
+
 #endif
diff --git a/include/boost/histogram/accumulators/weighted_sum.hpp b/include/boost/histogram/accumulators/weighted_sum.hpp
index 33ff8632..4754fdc0 100644
--- a/include/boost/histogram/accumulators/weighted_sum.hpp
+++ b/include/boost/histogram/accumulators/weighted_sum.hpp
@@ -7,7 +7,8 @@
 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_SUM_HPP
 #define BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_SUM_HPP
 
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
+#include <type_traits>
 
 namespace boost {
 namespace histogram {
@@ -88,4 +89,11 @@ private:
 } // namespace histogram
 } // namespace boost
 
+namespace std {
+template <class T, class U>
+struct common_type<boost::histogram::accumulators::weighted_sum<T>,
+                   boost::histogram::accumulators::weighted_sum<U>> {
+  using type = boost::histogram::accumulators::weighted_sum<common_type_t<T, U>>;
+};
+} // namespace std
 #endif
diff --git a/include/boost/histogram/adaptive_storage.hpp b/include/boost/histogram/adaptive_storage.hpp
index ece20135..82bbbd1b 100644
--- a/include/boost/histogram/adaptive_storage.hpp
+++ b/include/boost/histogram/adaptive_storage.hpp
@@ -15,7 +15,7 @@
 #include <boost/cstdint.hpp>
 #include <boost/histogram/detail/buffer.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/histogram/weight.hpp>
 #include <boost/iterator/iterator_facade.hpp>
 #include <boost/mp11.hpp>
diff --git a/include/boost/histogram/algorithm/project.hpp b/include/boost/histogram/algorithm/project.hpp
index 005fb0ad..61ee3bb0 100644
--- a/include/boost/histogram/algorithm/project.hpp
+++ b/include/boost/histogram/algorithm/project.hpp
@@ -23,24 +23,24 @@ namespace histogram {
 namespace algorithm {
 
 /**
-  Returns a lower-dimensional histogram, summing over removed axes.
+  Returns a lower-dimensional grid, summing over removed axes.
 
   Arguments are the source histogram and compile-time numbers, the remaining indices of
   the axes. Returns a new histogram which only contains the subset of axes. The source
   histogram is summed over the removed axes.
 */
-template <typename A, typename S, unsigned N, typename... Ns>
-auto project(const histogram<A, S>& hist, std::integral_constant<unsigned, N> n,
-             Ns... ns) {
+template <template <class, class> class Grid, class A, class S, unsigned N,
+          typename... Ns>
+auto project(const Grid<A, S>& grid, std::integral_constant<unsigned, N> n, Ns... ns) {
   using LN = mp11::mp_list<decltype(n), Ns...>;
   static_assert(mp11::mp_is_set<LN>::value, "indices must be unique");
 
-  auto axes = detail::make_sub_axes(unsafe_access::axes(hist), n, ns...);
-  auto result = histogram<decltype(axes), S>(
-      std::move(axes), detail::make_default(unsafe_access::storage(hist)));
+  auto axes = detail::make_sub_axes(unsafe_access::axes(grid), n, ns...);
+  auto result = Grid<decltype(axes), S>(
+      std::move(axes), detail::make_default(unsafe_access::storage(grid)));
 
-  detail::axes_buffer<typename decltype(result)::axes_type, int> idx(result.rank());
-  for (auto x : indexed(hist, true)) {
+  detail::axes_buffer<decltype(axes), int> idx(result.rank());
+  for (auto x : indexed(grid, true)) {
     auto i = idx.begin();
     mp11::mp_for_each<LN>([&i, &x](auto I) { *i++ = x[I]; });
     unsafe_access::add_value(result, idx, *x);
@@ -54,27 +54,27 @@ auto project(const histogram<A, S>& hist, std::integral_constant<unsigned, N> n,
   This version accepts a source histogram and an iterable range containing the remaining
   indices.
 */
-template <typename A, typename S, typename Iterable,
-          typename = detail::requires_iterable<Iterable>>
-auto project(const histogram<A, S>& hist, const Iterable& c) {
-  static_assert(detail::is_axis_vector<A>::value,
-                "dynamic version of project requires a histogram with dynamic axis");
+template <template <class, class> class Grid, class A, class S, class Iterable,
+          class = detail::requires_iterable<Iterable>>
+auto project(const Grid<A, S>& grid, const Iterable& c) {
+  static_assert(detail::is_sequence_of_any_axis<A>::value,
+                "dynamic version of project requires a grid with dynamic axis");
 
-  const auto& hist_axes = unsafe_access::axes(hist);
-  auto axes = detail::make_default(hist_axes);
+  const auto& old_axes = unsafe_access::axes(grid);
+  auto axes = detail::make_default(old_axes);
   axes.reserve(c.size());
-  detail::axes_buffer<A, bool> seen(hist_axes.size(), false);
+  detail::axes_buffer<A, bool> seen(old_axes.size(), false);
   for (auto d : c) {
     if (seen[d]) BOOST_THROW_EXCEPTION(std::invalid_argument("indices must be unique"));
     seen[d] = true;
-    axes.emplace_back(hist_axes[d]);
+    axes.emplace_back(old_axes[d]);
   }
 
-  auto result = histogram<A, S>(std::move(axes),
-                                detail::make_default(unsafe_access::storage(hist)));
+  auto result =
+      Grid<A, S>(std::move(axes), detail::make_default(unsafe_access::storage(grid)));
 
   detail::axes_buffer<A, int> idx(result.rank());
-  for (auto x : indexed(hist, true)) {
+  for (auto x : indexed(grid, true)) {
     auto i = idx.begin();
     for (auto d : c) *i++ = x[d];
     unsafe_access::add_value(result, idx, *x);
diff --git a/include/boost/histogram/algorithm/reduce.hpp b/include/boost/histogram/algorithm/reduce.hpp
index cd0539ab..774e9182 100644
--- a/include/boost/histogram/algorithm/reduce.hpp
+++ b/include/boost/histogram/algorithm/reduce.hpp
@@ -10,7 +10,7 @@
 #include <boost/assert.hpp>
 #include <boost/histogram/detail/axes.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/histogram/indexed.hpp>
 #include <boost/histogram/unsafe_access.hpp>
 
@@ -66,14 +66,16 @@ inline reduce_option_type shrink(double lower, double upper) {
 /// Convenience overload for single axis.
 inline reduce_option_type rebin(unsigned merge) { return rebin(0, merge); }
 
-template <typename A, typename S, typename C, typename = detail::requires_iterable<C>>
-histogram<A, S> reduce(const histogram<A, S>& hist, const C& options) {
+template <typename Grid, typename C, typename = detail::requires_iterable<C>>
+decltype(auto) reduce(const Grid& grid, const C& options) {
+  using axes_type = typename Grid::axes_type;
+
   struct option_item : reduce_option_type {
     int begin, end;
     bool is_set() const noexcept { return reduce_option_type::merge > 0; }
   };
 
-  auto options_internal = detail::axes_buffer<A, option_item>(hist.rank());
+  auto options_internal = detail::axes_buffer<axes_type, option_item>(grid.rank());
   for (const auto& o : options) {
     auto& oi = options_internal[o.iaxis];
     if (oi.is_set()) // did we already set the option for this axis?
@@ -83,10 +85,10 @@ histogram<A, S> reduce(const histogram<A, S>& hist, const C& options) {
     oi.merge = o.merge;
   }
 
-  auto axes = detail::make_empty_axes(unsafe_access::axes(hist));
+  auto axes = detail::make_empty_axes(unsafe_access::axes(grid));
 
   unsigned iaxis = 0;
-  hist.for_each_axis([&](const auto& a) {
+  grid.for_each_axis([&](const auto& a) {
     using T = detail::unqual<decltype(a)>;
 
     auto& o = options_internal[iaxis];
@@ -111,11 +113,10 @@ histogram<A, S> reduce(const histogram<A, S>& hist, const C& options) {
     ++iaxis;
   });
 
-  auto result = histogram<A, S>(std::move(axes),
-                                detail::make_default(unsafe_access::storage(hist)));
+  auto result = Grid(std::move(axes), detail::make_default(unsafe_access::storage(grid)));
 
-  detail::axes_buffer<A, int> idx(hist.rank());
-  for (auto x : indexed(hist, true)) {
+  detail::axes_buffer<axes_type, int> idx(grid.rank());
+  for (auto x : indexed(grid, true)) {
     auto i = idx.begin();
     auto o = options_internal.begin();
     for (auto j : x) {
@@ -136,11 +137,10 @@ histogram<A, S> reduce(const histogram<A, S>& hist, const C& options) {
   return result;
 }
 
-template <typename A, typename S, typename... Ts>
-histogram<A, S> reduce(const histogram<A, S>& h, const reduce_option_type& t,
-                       Ts&&... ts) {
+template <class Grid, class... Ts>
+decltype(auto) reduce(const Grid& grid, const reduce_option_type& t, Ts&&... ts) {
   // this must be in one line, because any of the ts could be a temporary
-  return reduce(h, std::initializer_list<reduce_option_type>{t, ts...});
+  return reduce(grid, std::initializer_list<reduce_option_type>{t, ts...});
 }
 
 } // namespace algorithm
diff --git a/include/boost/histogram/algorithm/sum.hpp b/include/boost/histogram/algorithm/sum.hpp
index e290d454..fa04f5d5 100644
--- a/include/boost/histogram/algorithm/sum.hpp
+++ b/include/boost/histogram/algorithm/sum.hpp
@@ -8,21 +8,22 @@
 #define BOOST_HISTOGRAM_ALGORITHM_SUM_HPP
 
 #include <boost/histogram/accumulators/sum.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <numeric>
 #include <type_traits>
 
 namespace boost {
 namespace histogram {
 namespace algorithm {
-template <typename A, typename S,
-          typename ReturnType = std::conditional_t<
-              std::is_arithmetic<typename histogram<A, S>::value_type>::value, double,
-              typename histogram<A, S>::value_type>,
-          typename InternalSum = std::conditional_t<
-              std::is_arithmetic<typename histogram<A, S>::value_type>::value,
-              accumulators::sum<double>, typename histogram<A, S>::value_type>>
-ReturnType sum(const histogram<A, S>& h) {
+template <
+    typename A, typename S,
+    typename ReturnType =
+        std::conditional_t<std::is_arithmetic<typename grid<A, S>::value_type>::value,
+                           double, typename grid<A, S>::value_type>,
+    typename InternalSum =
+        std::conditional_t<std::is_arithmetic<typename grid<A, S>::value_type>::value,
+                           accumulators::sum<double>, typename grid<A, S>::value_type>>
+ReturnType sum(const grid<A, S>& h) {
   InternalSum sum;
   for (auto x : h) sum += x;
   return sum;
diff --git a/include/boost/histogram/arithmetic_operators.hpp b/include/boost/histogram/arithmetic_operators.hpp
deleted file mode 100644
index 9c1582f7..00000000
--- a/include/boost/histogram/arithmetic_operators.hpp
+++ /dev/null
@@ -1,83 +0,0 @@
-// Copyright 2015-2017 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)
-
-#ifndef BOOST_HISTOGRAM_HISTOGRAM_ARITHMETIC_OPERATORS_HPP
-#define BOOST_HISTOGRAM_HISTOGRAM_ARITHMETIC_OPERATORS_HPP
-
-#include <boost/histogram/histogram_fwd.hpp>
-
-namespace boost {
-namespace histogram {
-
-template <typename A, typename S>
-histogram<A, S>&& operator+(histogram<A, S>&& a, const histogram<A, S>& b) {
-  a += b;
-  return std::move(a);
-}
-
-template <typename A, typename S>
-histogram<A, S>&& operator+(histogram<A, S>&& a, histogram<A, S>&& b) {
-  a += b;
-  return std::move(a);
-}
-
-template <typename A, typename S>
-histogram<A, S>&& operator+(const histogram<A, S>& a, histogram<A, S>&& b) {
-  b += a;
-  return std::move(b);
-}
-
-template <typename A, typename S>
-histogram<A, S> operator+(const histogram<A, S>& a,
-                          const histogram<A, S>& b) {
-  histogram<A, S> r(a);
-  r += b;
-  return r;
-}
-
-template <typename A, typename S>
-histogram<A, S>&& operator*(histogram<A, S>&& a, const double x) {
-  a *= x;
-  return std::move(a);
-}
-
-template <typename A, typename S>
-histogram<A, S>&& operator*(const double x, histogram<A, S>&& b) {
-  b *= x;
-  return std::move(b);
-}
-
-template <typename A, typename S>
-histogram<A, S> operator*(const histogram<A, S>& a, const double x) {
-  auto r = a;
-  r *= x;
-  return r;
-}
-
-template <typename A, typename S>
-histogram<A, S> operator*(const double x, const histogram<A, S>& b) {
-  auto r = b;
-  r *= x;
-  return r;
-}
-
-template <typename A, typename S>
-histogram<A, S>&& operator/(histogram<A, S>&& a, const double x) {
-  a /= x;
-  return std::move(a);
-}
-
-template <typename A, typename S>
-histogram<A, S> operator/(const histogram<A, S>& a, const double x) {
-  auto r = a;
-  r /= x;
-  return r;
-}
-
-} // namespace histogram
-} // namespace boost
-
-#endif
diff --git a/include/boost/histogram/axis/base.hpp b/include/boost/histogram/axis/base.hpp
index d0aeb774..20b181b0 100644
--- a/include/boost/histogram/axis/base.hpp
+++ b/include/boost/histogram/axis/base.hpp
@@ -10,7 +10,7 @@
 #include <boost/histogram/detail/cat.hpp>
 #include <boost/histogram/detail/compressed_pair.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/throw_exception.hpp>
 #include <limits>
 #include <stdexcept>
diff --git a/include/boost/histogram/axis/category.hpp b/include/boost/histogram/axis/category.hpp
index 8884e540..6e6983f0 100644
--- a/include/boost/histogram/axis/category.hpp
+++ b/include/boost/histogram/axis/category.hpp
@@ -15,7 +15,7 @@
 #include <boost/histogram/detail/buffer.hpp>
 #include <boost/histogram/detail/compressed_pair.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/throw_exception.hpp>
 #include <memory>
 #include <stdexcept>
diff --git a/include/boost/histogram/axis/integer.hpp b/include/boost/histogram/axis/integer.hpp
index 7163a8ad..a163f432 100644
--- a/include/boost/histogram/axis/integer.hpp
+++ b/include/boost/histogram/axis/integer.hpp
@@ -11,7 +11,7 @@
 #include <boost/histogram/axis/base.hpp>
 #include <boost/histogram/axis/iterator.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/throw_exception.hpp>
 #include <cmath>
 #include <limits>
diff --git a/include/boost/histogram/axis/ostream_operators.hpp b/include/boost/histogram/axis/ostream_operators.hpp
index 7e31fc4c..01affae5 100644
--- a/include/boost/histogram/axis/ostream_operators.hpp
+++ b/include/boost/histogram/axis/ostream_operators.hpp
@@ -16,6 +16,7 @@
 #include <boost/throw_exception.hpp>
 #include <iomanip>
 #include <ostream>
+#include <stdexcept>
 #include <type_traits>
 
 namespace boost {
diff --git a/include/boost/histogram/axis/regular.hpp b/include/boost/histogram/axis/regular.hpp
index 674c3e4a..d3514be3 100644
--- a/include/boost/histogram/axis/regular.hpp
+++ b/include/boost/histogram/axis/regular.hpp
@@ -12,7 +12,7 @@
 #include <boost/histogram/axis/interval_view.hpp>
 #include <boost/histogram/axis/iterator.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/mp11.hpp>
 #include <boost/throw_exception.hpp>
 #include <cmath>
diff --git a/include/boost/histogram/axis/traits.hpp b/include/boost/histogram/axis/traits.hpp
index f12f3ceb..965ba5e4 100644
--- a/include/boost/histogram/axis/traits.hpp
+++ b/include/boost/histogram/axis/traits.hpp
@@ -10,8 +10,9 @@
 #include <boost/core/typeinfo.hpp>
 #include <boost/histogram/detail/cat.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/throw_exception.hpp>
+#include <stdexcept>
 #include <utility>
 
 namespace boost {
diff --git a/include/boost/histogram/axis/variable.hpp b/include/boost/histogram/axis/variable.hpp
index 63410833..8218057e 100644
--- a/include/boost/histogram/axis/variable.hpp
+++ b/include/boost/histogram/axis/variable.hpp
@@ -16,7 +16,7 @@
 #include <boost/histogram/detail/buffer.hpp>
 #include <boost/histogram/detail/compressed_pair.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/throw_exception.hpp>
 #include <cmath>
 #include <limits>
diff --git a/include/boost/histogram/axis/variant.hpp b/include/boost/histogram/axis/variant.hpp
index 770f91a1..1b77b4fc 100644
--- a/include/boost/histogram/axis/variant.hpp
+++ b/include/boost/histogram/axis/variant.hpp
@@ -13,7 +13,7 @@
 #include <boost/histogram/axis/traits.hpp>
 #include <boost/histogram/detail/cat.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/mp11.hpp>
 #include <boost/throw_exception.hpp>
 #include <boost/variant.hpp>
diff --git a/include/boost/histogram/detail/axes.hpp b/include/boost/histogram/detail/axes.hpp
index 6863603e..6d0a145b 100644
--- a/include/boost/histogram/detail/axes.hpp
+++ b/include/boost/histogram/detail/axes.hpp
@@ -13,7 +13,7 @@
 #include <boost/histogram/axis/traits.hpp>
 #include <boost/histogram/axis/variant.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/mp11.hpp>
 #include <boost/throw_exception.hpp>
 #include <stdexcept>
diff --git a/include/boost/histogram/detail/common_type.hpp b/include/boost/histogram/detail/common_type.hpp
new file mode 100644
index 00000000..93ce33a1
--- /dev/null
+++ b/include/boost/histogram/detail/common_type.hpp
@@ -0,0 +1,62 @@
+// Copyright 2015-2018 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)
+
+#ifndef BOOST_HISTOGRAM_DETAIL_COMMON_TYPE_HPP
+#define BOOST_HISTOGRAM_DETAIL_COMMON_TYPE_HPP
+
+#include <boost/histogram/detail/meta.hpp>
+#include <boost/histogram/fwd.hpp>
+#include <boost/mp11.hpp>
+#include <tuple>
+#include <type_traits>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+template <class T, class U>
+// clang-format off
+using common_axes = mp11::mp_cond<
+  std::is_same<T, U>, T,
+  is_tuple<T>, T,
+  is_tuple<U>, U,
+  is_sequence_of_axis<T>, T,
+  is_sequence_of_axis<U>, U,
+  mp11::mp_true, T
+>;
+// clang-format on
+
+template <class T, class U>
+struct common_storage_impl;
+
+template <class T, class U>
+struct common_storage_impl<storage_adaptor<T>, storage_adaptor<U>> {
+  using V = std::common_type_t<mp11::mp_front<T>, mp11::mp_front<U>>;
+  using type = storage_adaptor<mp11::mp_replace_front<T, V>>;
+};
+
+template <class T, class A>
+struct common_storage_impl<storage_adaptor<T>, adaptive_storage<A>> {
+  using U = typename adaptive_storage<A>::value_type;
+  using V = std::common_type_t<mp11::mp_front<T>, U>;
+  using type = storage_adaptor<mp11::mp_replace_front<T, V>>;
+};
+
+template <class T, class A>
+struct common_storage_impl<adaptive_storage<A>, storage_adaptor<T>>
+    : common_storage_impl<storage_adaptor<T>, adaptive_storage<A>> {};
+
+template <class A1, class A2>
+struct common_storage_impl<adaptive_storage<A1>, adaptive_storage<A2>> {
+  using type = adaptive_storage<A1>;
+};
+
+template <class A, class B>
+using common_storage = typename common_storage_impl<A, B>::type;
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif
diff --git a/include/boost/histogram/detail/is_set.hpp b/include/boost/histogram/detail/is_set.hpp
deleted file mode 100644
index b87eafd5..00000000
--- a/include/boost/histogram/detail/is_set.hpp
+++ /dev/null
@@ -1,29 +0,0 @@
-// Copyright 2018 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)
-
-#ifndef BOOST_HISTOGRAM_DETAIL_IS_SET_HPP
-#define BOOST_HISTOGRAM_DETAIL_IS_SET_HPP
-
-#include <boost/container/flat_set.hpp>
-#include <boost/container/static_vector.hpp>
-#include <boost/histogram/detail/meta.hpp>
-
-namespace boost {
-namespace histogram {
-namespace detail {
-template <class Iterator>
-// TODO: remove this
-bool is_set(Iterator begin, Iterator end) {
-  using T = iterator_value_type<Iterator>;
-  using C = boost::container::static_vector<T, BOOST_HISTOGRAM_DETAIL_AXES_LIMIT>;
-  boost::container::flat_set<T, std::less<T>, C> s(begin, end);
-  return static_cast<std::size_t>(std::distance(begin, end)) == s.size();
-}
-} // namespace detail
-} // namespace histogram
-} // namespace boost
-
-#endif
diff --git a/include/boost/histogram/detail/meta.hpp b/include/boost/histogram/detail/meta.hpp
index 545a1426..13021293 100644
--- a/include/boost/histogram/detail/meta.hpp
+++ b/include/boost/histogram/detail/meta.hpp
@@ -17,7 +17,7 @@
 #if BOOST_WORKAROUND(BOOST_GCC, >= 60000)
 #pragma GCC diagnostic pop
 #endif
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/mp11.hpp>
 #include <functional>
 #include <iterator>
@@ -58,12 +58,6 @@ using mp_set_union = mp11::mp_apply_q<mp11::mp_bind_front<mp11::mp_set_push_back
 template <typename L>
 using mp_last = mp11::mp_at_c<L, (mp_size<L>::value - 1)>;
 
-template <typename T>
-using container_value_type = mp11::mp_first<unqual<T>>;
-
-template <typename T>
-using iterator_value_type = typename std::iterator_traits<T>::value_type;
-
 template <typename T>
 using args_type = mp11::mp_if<std::is_member_function_pointer<T>,
                               mp11::mp_pop_front<boost::callable_traits::args_t<T>>,
@@ -241,14 +235,18 @@ template <typename T>
 using is_axis_variant = typename is_axis_variant_impl<T>::type;
 
 template <typename T>
-using is_axis_or_axis_variant = mp11::mp_or<is_axis<T>, is_axis_variant<T>>;
+using is_any_axis = mp11::mp_or<is_axis<T>, is_axis_variant<T>>;
 
 template <typename T>
-using is_axes = is_axis_or_axis_variant<mp11::mp_first<T>>;
+using is_sequence_of_axis = mp11::mp_and<is_iterable<T>, is_axis<mp11::mp_first<T>>>;
 
-template <typename T, typename U = container_value_type<T>>
-using is_axis_vector =
-    mp11::mp_all<is_indexable_container<unqual<T>>, is_axis_or_axis_variant<U>>;
+template <typename T>
+using is_sequence_of_axis_variant =
+    mp11::mp_and<is_iterable<T>, is_axis_variant<mp11::mp_first<T>>>;
+
+template <typename T>
+using is_sequence_of_any_axis =
+    mp11::mp_and<is_iterable<T>, is_any_axis<mp11::mp_first<T>>>;
 
 template <typename T>
 struct is_weight_impl : std::false_type {};
@@ -278,13 +276,21 @@ struct requires_iterable {};
 template <typename T, typename = mp11::mp_if<is_axis<unqual<T>>, void>>
 struct requires_axis {};
 
-template <typename T, typename = mp11::mp_if<is_axis_or_axis_variant<unqual<T>>, void>>
-struct requires_axis_or_axis_variant {};
+template <typename T, typename = mp11::mp_if<is_any_axis<unqual<T>>, void>>
+struct requires_any_axis {};
 
-template <typename T, typename = mp11::mp_if<is_axis_vector<unqual<T>>, void>>
-struct requires_axis_vector {};
+template <typename T, typename = mp11::mp_if<is_sequence_of_axis<unqual<T>>, void>>
+struct requires_sequence_of_axis {};
 
-template <typename T, typename = mp11::mp_if<is_axes<unqual<T>>, void>>
+template <typename T,
+          typename = mp11::mp_if<is_sequence_of_axis_variant<unqual<T>>, void>>
+struct requires_sequence_of_axis_variant {};
+
+template <typename T, typename = mp11::mp_if<is_sequence_of_any_axis<unqual<T>>, void>>
+struct requires_sequence_of_any_axis {};
+
+template <typename T,
+          typename = mp11::mp_if<is_any_axis<mp11::mp_first<unqual<T>>>, void>>
 struct requires_axes {};
 
 template <typename T>
diff --git a/include/boost/histogram/histogram_fwd.hpp b/include/boost/histogram/fwd.hpp
similarity index 93%
rename from include/boost/histogram/histogram_fwd.hpp
rename to include/boost/histogram/fwd.hpp
index 97438f9f..3814cb59 100644
--- a/include/boost/histogram/histogram_fwd.hpp
+++ b/include/boost/histogram/fwd.hpp
@@ -4,11 +4,11 @@
 // (See accompanying file LICENSE_1_0.txt
 // or copy at http://www.boost.org/LICENSE_1_0.txt)
 
-#ifndef BOOST_HISTOGRAM_HISTOGRAM_FWD_HPP
-#define BOOST_HISTOGRAM_HISTOGRAM_FWD_HPP
+#ifndef BOOST_fwd_HPP
+#define BOOST_fwd_HPP
 
-#include <boost/config.hpp>
 #include <boost/container/container_fwd.hpp> // for string and new_allocator
+#include <boost/histogram/attribute.hpp>     // for BOOST_HISTOGRAM_NODISCARD
 
 // Why boost containers as defaults and not std containers?
 // - std::vector does not work with incomplete types, but boost::container::vector does
@@ -109,8 +109,11 @@ using profile_storage = storage_adaptor<boost::container::vector<accumulators::m
 using weighted_profile_storage =
     storage_adaptor<boost::container::vector<accumulators::weighted_mean<>>>;
 
+template <class Axes, class Storage>
+class BOOST_HISTOGRAM_NODISCARD grid;
+
 template <class Axes, class Storage = default_storage>
-class histogram;
+class BOOST_HISTOGRAM_NODISCARD histogram;
 } // namespace histogram
 } // namespace boost
 
diff --git a/include/boost/histogram/grid.hpp b/include/boost/histogram/grid.hpp
new file mode 100644
index 00000000..3fe3d711
--- /dev/null
+++ b/include/boost/histogram/grid.hpp
@@ -0,0 +1,253 @@
+// Copyright 2015-2018 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)
+
+#ifndef BOOST_HISTOGRAM_GRID_HPP
+#define BOOST_HISTOGRAM_GRID_HPP
+
+#include <boost/histogram/detail/axes.hpp>
+#include <boost/histogram/detail/meta.hpp>
+#include <boost/histogram/fwd.hpp>
+#include <boost/mp11.hpp>
+#include <boost/throw_exception.hpp>
+#include <stdexcept>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+namespace boost {
+namespace histogram {
+
+template <class Axes, class Storage>
+class grid {
+  static_assert(mp11::mp_size<Axes>::value > 0, "at least one axis required");
+
+public:
+  using axes_type = Axes;
+  using storage_type = Storage;
+  using value_type = typename storage_type::value_type;
+  // typedefs for boost::range_iterator
+  using iterator = decltype(std::declval<storage_type&>().begin());
+  using const_iterator = decltype(std::declval<const storage_type&>().begin());
+
+  grid() = default;
+  grid(const grid& rhs) = default;
+  grid(grid&& rhs) = default;
+  grid& operator=(const grid& rhs) = default;
+  grid& operator=(grid&& rhs) = default;
+
+  template <class A, class S>
+  explicit grid(const grid<A, S>& rhs) : storage_(rhs.storage_) {
+    detail::axes_assign(axes_, rhs.axes_);
+  }
+
+  template <class A, class S>
+  grid& operator=(const grid<A, S>& rhs) {
+    storage_ = rhs.storage_;
+    detail::axes_assign(axes_, rhs.axes_);
+    return *this;
+  }
+
+  template <class A, class S>
+  explicit grid(A&& a, S&& s) : axes_(std::forward<A>(a)), storage_(std::forward<S>(s)) {
+    storage_.reset(detail::bincount(axes_));
+  }
+
+  template <class A, class S>
+  bool operator==(const grid<A, S>& rhs) const noexcept {
+    return detail::axes_equal(axes_, rhs.axes_) && storage_ == rhs.storage_;
+  }
+
+  template <class A, class S>
+  bool operator!=(const grid<A, S>& rhs) const noexcept {
+    return !operator==(rhs);
+  }
+
+  template <class A, class S>
+  grid& operator+=(const grid<A, S>& rhs) {
+    if (!detail::axes_equal(axes_, rhs.axes_))
+      BOOST_THROW_EXCEPTION(std::invalid_argument("axes of histograms differ"));
+    storage_ += rhs.storage_;
+    return *this;
+  }
+
+  grid& operator*=(const double x) {
+    storage_ *= x;
+    return *this;
+  }
+
+  grid& operator/=(const double x) {
+    storage_ *= 1.0 / x;
+    return *this;
+  }
+
+  /// Number of axes (dimensions)
+  unsigned rank() const noexcept { return detail::axes_size(axes_); }
+
+  /// Total number of bins (including underflow/overflow)
+  std::size_t size() const noexcept { return storage_.size(); }
+
+  /// Reset values to zero
+  void reset() { storage_.reset(storage_.size()); }
+
+  /// Get N-th axis
+  template <unsigned N>
+  decltype(auto) axis(std::integral_constant<unsigned, N>) {
+    detail::rank_check(axes_, N);
+    return detail::axis_get<N>(axes_);
+  }
+
+  /// Get N-th axis (const version)
+  template <unsigned N>
+  decltype(auto) axis(std::integral_constant<unsigned, N>) const {
+    detail::rank_check(axes_, N);
+    return detail::axis_get<N>(axes_);
+  }
+
+  /// Get first axis (convenience for 1-d histograms)
+  decltype(auto) axis() { return axis(std::integral_constant<unsigned, 0>()); }
+
+  /// Get first axis (convenience for 1-d histograms, const version)
+  decltype(auto) axis() const { return axis(std::integral_constant<unsigned, 0>()); }
+
+  /// Get N-th axis with runtime index
+  decltype(auto) axis(unsigned i) {
+    detail::rank_check(axes_, i);
+    return detail::axis_get(axes_, i);
+  }
+
+  /// Get N-th axis with runtime index (const version)
+  decltype(auto) axis(unsigned i) const {
+    detail::rank_check(axes_, i);
+    return detail::axis_get(axes_, i);
+  }
+
+  /// Apply unary functor/function to each axis
+  template <class Unary>
+  void for_each_axis(Unary&& unary) const {
+    detail::for_each_axis(axes_, std::forward<Unary>(unary));
+  }
+
+  /// Access value at indices
+  template <typename... Ts>
+  decltype(auto) at(const Ts&... ts) {
+    return at(std::forward_as_tuple(ts...));
+  }
+
+  /// Access value at indices (const version)
+  template <typename... Ts>
+  decltype(auto) at(const Ts&... ts) const {
+    return at(std::forward_as_tuple(ts...));
+  }
+
+  /// Access value at index tuple
+  template <typename... Ts>
+  decltype(auto) at(const std::tuple<Ts...>& t) {
+    const auto idx = detail::at_impl(axes_, t);
+    if (!idx) BOOST_THROW_EXCEPTION(std::out_of_range("indices out of bounds"));
+    return storage_[*idx];
+  }
+
+  /// Access value at index tuple (const version)
+  template <typename... Ts>
+  decltype(auto) at(const std::tuple<Ts...>& t) const {
+    const auto idx = detail::at_impl(axes_, t);
+    if (!idx) BOOST_THROW_EXCEPTION(std::out_of_range("indices out of bounds"));
+    return storage_[*idx];
+  }
+
+  /// Access value at index iterable
+  template <class Iterable, class = detail::requires_iterable<Iterable>>
+  decltype(auto) at(const Iterable& c) {
+    const auto idx = detail::at_impl(axes_, c);
+    if (!idx) BOOST_THROW_EXCEPTION(std::out_of_range("indices out of bounds"));
+    return storage_[*idx];
+  }
+
+  /// Access value at index iterable (const version)
+  template <class Iterable, class = detail::requires_iterable<Iterable>>
+  decltype(auto) at(const Iterable& c) const {
+    const auto idx = detail::at_impl(axes_, c);
+    if (!idx) BOOST_THROW_EXCEPTION(std::out_of_range("indices out of bounds"));
+    return storage_[*idx];
+  }
+
+  /// Access value at index (number for rank=1 or index tuple|iterable)
+  template <class T>
+  decltype(auto) operator[](const T& t) {
+    return at(t);
+  }
+
+  /// Access value at index (number for rank=1 or index tuple|iterable, const version)
+  template <class T>
+  decltype(auto) operator[](const T& t) const {
+    return at(t);
+  }
+
+  // iterator begin() noexcept { return storage_.begin(); }
+  // iterator end() noexcept { return storage_.end(); }
+
+  const_iterator begin() const noexcept { return storage_.begin(); }
+  const_iterator end() const noexcept { return storage_.end(); }
+
+protected:
+  axes_type axes_;
+  storage_type storage_;
+
+  template <class A, class S>
+  friend class grid;
+  friend struct unsafe_access;
+};
+
+template <class A, class S>
+auto operator*(grid<A, S>&& g, double x) {
+  return g *= x;
+}
+
+template <class A, class S>
+auto operator*(double x, grid<A, S>&& g) {
+  return g *= x;
+}
+
+template <class A, class S>
+auto operator*(const grid<A, S>& g, double x) {
+  auto copy = g;
+  return copy *= x;
+}
+
+template <class A, class S>
+auto operator*(double x, const grid<A, S>& g) {
+  auto copy = g;
+  return copy *= x;
+}
+
+template <class A, class S>
+auto operator/(grid<A, S>&& g, double x) {
+  return g *= 1.0 / x;
+}
+
+template <class A, class S>
+auto operator/(const grid<A, S>& g, double x) {
+  auto copy = g;
+  return copy *= 1.0 / x;
+}
+
+template <class A1, class S1, class A2, class S2>
+auto operator+(const grid<A1, S1>& a, const grid<A2, S2>& b) {
+  auto r = grid<detail::common_axes<A1, A2>, detail::common_storage<S1, S2>>(a);
+  return r += b;
+}
+
+#if __cpp_deduction_guides >= 201606
+
+template <class Axes, class Storage>
+grid(Axes&& axes, Storage&& storage)->grid<detail::unqual<Axes>, detail::unqual<Storage>>;
+
+#endif
+
+} // namespace histogram
+} // namespace boost
+
+#endif
diff --git a/include/boost/histogram/histogram.hpp b/include/boost/histogram/histogram.hpp
index 0361936c..3c55683d 100644
--- a/include/boost/histogram/histogram.hpp
+++ b/include/boost/histogram/histogram.hpp
@@ -7,191 +7,114 @@
 #ifndef BOOST_HISTOGRAM_HISTOGRAM_HPP
 #define BOOST_HISTOGRAM_HISTOGRAM_HPP
 
-#include <algorithm>
-#include <boost/assert.hpp>
-#include <boost/histogram/arithmetic_operators.hpp>
-#include <boost/histogram/attribute.hpp>
 #include <boost/histogram/detail/axes.hpp>
+#include <boost/histogram/detail/common_type.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
-#include <boost/mp11.hpp>
-#include <boost/throw_exception.hpp>
-#include <functional>
+#include <boost/histogram/fwd.hpp>
+#include <boost/histogram/grid.hpp>
 #include <tuple>
 #include <type_traits>
-#include <utility>
 
 namespace boost {
 namespace histogram {
 
-template <typename Axes, typename Storage>
-class BOOST_HISTOGRAM_NODISCARD histogram {
-  static_assert(mp11::mp_size<Axes>::value > 0, "at least one axis required");
+template <class Axes, class Storage>
+class histogram : public grid<Axes, Storage> {
+  using grid_type = grid<Axes, Storage>;
 
 public:
-  using axes_type = Axes;
-  using storage_type = Storage;
-  using value_type = typename storage_type::value_type;
-  // typedefs for boost::range_iterator
-  using iterator = decltype(std::declval<storage_type&>().begin());
-  using const_iterator = decltype(std::declval<const storage_type&>().begin());
-
   histogram() = default;
   histogram(const histogram& rhs) = default;
   histogram(histogram&& rhs) = default;
   histogram& operator=(const histogram& rhs) = default;
   histogram& operator=(histogram&& rhs) = default;
 
-  template <typename A, typename S>
-  explicit histogram(const histogram<A, S>& rhs) : storage_(rhs.storage_) {
-    detail::axes_assign(axes_, rhs.axes_);
+  template <class A, class S>
+  explicit histogram(const histogram<A, S>& rhs) : grid_type(rhs) {}
+
+  template <class A, class = detail::requires_axes<A>>
+  explicit histogram(A&& a)
+      : grid_type(std::forward<A>(a), typename grid_type::storage_type()) {}
+
+  template <class A, class S>
+  histogram(A&& a, S&& s) : grid_type(std::forward<A>(a), std::forward<S>(s)) {}
+
+  template <class A, class S>
+  histogram& operator=(histogram<A, S>&& rhs) {
+    grid_type::operator=(std::move(rhs));
+    return *this;
   }
 
-  template <typename A, typename S>
+  template <class A, class S>
   histogram& operator=(const histogram<A, S>& rhs) {
-    storage_ = rhs.storage_;
-    detail::axes_assign(axes_, rhs.axes_);
+    grid_type::operator=(rhs);
     return *this;
   }
 
-  template <typename A, typename = detail::requires_axes<A>>
-  explicit histogram(A&& a) : axes_(std::forward<A>(a)) {
-    storage_.reset(detail::bincount(axes_));
-  }
-
-  template <typename A, typename S>
-  explicit histogram(A&& a, S&& s)
-      : axes_(std::forward<A>(a)), storage_(std::forward<S>(s)) {
-    storage_.reset(detail::bincount(axes_));
-  }
-
-  template <typename A, typename S>
-  bool operator==(const histogram<A, S>& rhs) const noexcept {
-    return detail::axes_equal(axes_, rhs.axes_) && storage_ == rhs.storage_;
-  }
-
-  template <typename A, typename S>
-  bool operator!=(const histogram<A, S>& rhs) const noexcept {
-    return !operator==(rhs);
-  }
-
-  template <typename A, typename S>
+  /// Add values of another histogram
+  template <class A, class S>
   histogram& operator+=(const histogram<A, S>& rhs) {
-    if (!detail::axes_equal(axes_, rhs.axes_))
-      BOOST_THROW_EXCEPTION(std::invalid_argument("axes of histograms differ"));
-    storage_ += rhs.storage_;
+    grid_type::operator+=(rhs);
     return *this;
   }
 
-  histogram& operator*=(const double x) {
-    storage_ *= x;
-    return *this;
-  }
-
-  histogram& operator/=(const double x) {
-    storage_ *= 1.0 / x;
-    return *this;
-  }
-
-  /// Number of axes (dimensions)
-  unsigned rank() const noexcept { return detail::axes_size(axes_); }
-
-  /// Total number of bins (including underflow/overflow)
-  std::size_t size() const noexcept { return storage_.size(); }
-
-  /// Reset values to zero
-  void reset() { storage_.reset(storage_.size()); }
-
-  /// Get N-th axis (const version)
-  template <unsigned N>
-  decltype(auto) axis(std::integral_constant<unsigned, N>) const {
-    detail::rank_check(axes_, N);
-    return detail::axis_get<N>(axes_);
-  }
-
-  /// Get N-th axis
-  template <unsigned N>
-  decltype(auto) axis(std::integral_constant<unsigned, N>) {
-    detail::rank_check(axes_, N);
-    return detail::axis_get<N>(axes_);
-  }
-
-  /// Get first axis (convenience for 1-d histograms, const version)
-  decltype(auto) axis() const { return axis(std::integral_constant<unsigned, 0>()); }
-
-  /// Get first axis (convenience for 1-d histograms)
-  decltype(auto) axis() { return axis(std::integral_constant<unsigned, 0>()); }
-
-  /// Get N-th axis with runtime index (const version)
-  decltype(auto) axis(unsigned i) const {
-    detail::rank_check(axes_, i);
-    return detail::axis_get(axes_, i);
-  }
-
-  /// Get N-th axis with runtime index
-  decltype(auto) axis(unsigned i) {
-    detail::rank_check(axes_, i);
-    return detail::axis_get(axes_, i);
-  }
-
-  /// Apply unary functor/function to each axis
-  template <typename Unary>
-  void for_each_axis(Unary&& unary) const {
-    detail::for_each_axis(axes_, std::forward<Unary>(unary));
-  }
-
   /// Fill histogram with values and optional weight or sample
-  template <typename... Ts>
+  template <class... Ts>
   void operator()(const Ts&... ts) {
     operator()(std::forward_as_tuple(ts...));
   }
 
   /// Fill histogram with value tuple and optional weight or sample
-  template <typename... Ts>
+  template <class... Ts>
   void operator()(const std::tuple<Ts...>& t) {
-    detail::fill_impl(storage_, axes_, t);
+    detail::fill_impl(grid_type::storage_, grid_type::axes_, t);
   }
 
-  /// Access value at indices
+  /// Access value at indices (const version)
   template <typename... Ts>
   decltype(auto) at(const Ts&... ts) const {
     return at(std::forward_as_tuple(ts...));
   }
 
-  /// Access value at index tuple
+  /// Access value at index tuple (const version)
   template <typename... Ts>
   decltype(auto) at(const std::tuple<Ts...>& t) const {
-    const auto idx = detail::at_impl(axes_, t);
-    if (!idx) BOOST_THROW_EXCEPTION(std::out_of_range("indices out of bounds"));
-    return storage_[*idx];
+    return grid_type::at(t);
   }
 
-  /// Access value at index iterable
-  template <typename Iterable, typename = detail::requires_iterable<Iterable>>
+  /// Access value at index iterable (const version)
+  template <class Iterable, class = detail::requires_iterable<Iterable>>
   decltype(auto) at(const Iterable& c) const {
-    const auto idx = detail::at_impl(axes_, c);
-    if (!idx) BOOST_THROW_EXCEPTION(std::out_of_range("indices out of bounds"));
-    return storage_[*idx];
+    return grid_type::at(c);
   }
 
-  /// Access value at index (number for rank=1 or index tuple|iterable)
-  template <typename T>
+  /// Access value at index (number for rank=1 or index tuple|iterable, const version)
+  template <class T>
   decltype(auto) operator[](const T& t) const {
     return at(t);
   }
 
-  const_iterator begin() const noexcept { return storage_.begin(); }
-  const_iterator end() const noexcept { return storage_.end(); }
-
-private:
-  axes_type axes_;
-  storage_type storage_;
-
-  template <typename A, typename S>
-  friend class histogram;
-  friend struct unsafe_access;
+  auto begin() const noexcept { return grid_type::begin(); }
+  auto end() const noexcept { return grid_type::end(); }
 };
 
+template <class A1, class S1, class A2, class S2>
+auto operator+(const histogram<A1, S1>& a, const histogram<A2, S2>& b) {
+  histogram<detail::common_axes<A1, A2>, detail::common_storage<S1, S2>> result = a;
+  result += b;
+  return result;
+}
+
+template <class A, class S>
+auto operator+(histogram<A, S>&& a, const histogram<A, S>& b) {
+  return a += b;
+}
+
+template <class A, class S>
+auto operator+(const histogram<A, S>& a, histogram<A, S>&& b) {
+  return b += a;
+}
+
 #if __cpp_deduction_guides >= 201606
 
 template <class Axes>
diff --git a/include/boost/histogram/indexed.hpp b/include/boost/histogram/indexed.hpp
index 2aa24f9f..7304280d 100644
--- a/include/boost/histogram/indexed.hpp
+++ b/include/boost/histogram/indexed.hpp
@@ -11,7 +11,7 @@
 #include <boost/histogram/axis/traits.hpp>
 #include <boost/histogram/detail/axes.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/iterator/iterator_adaptor.hpp>
 #include <boost/mp11.hpp>
 #include <type_traits>
diff --git a/include/boost/histogram/make_histogram.hpp b/include/boost/histogram/make_histogram.hpp
index 18994803..eb494506 100644
--- a/include/boost/histogram/make_histogram.hpp
+++ b/include/boost/histogram/make_histogram.hpp
@@ -46,22 +46,25 @@ auto make_weighted_histogram(T&& axis0, Ts&&... axis) {
 
 /// histogram factory from vector-like with custom storage
 template <typename StorageOrContainer, typename Iterable,
-          typename = detail::requires_axis_vector<Iterable>>
+          typename = detail::requires_sequence_of_any_axis<Iterable>>
 auto make_histogram_with(StorageOrContainer&& s, Iterable&& c) {
   using U = detail::unqual<StorageOrContainer>;
   using S = mp11::mp_if<detail::is_storage<U>, U, storage_adaptor<U>>;
-  return histogram<detail::unqual<Iterable>, S>(std::forward<Iterable>(c),
-                                                S(std::forward<StorageOrContainer>(s)));
+  using It = detail::unqual<Iterable>;
+  using A = mp11::mp_if<detail::is_indexable_container<It>, It,
+                        boost::container::vector<mp11::mp_first<It>>>;
+  return histogram<A, S>(std::forward<Iterable>(c),
+                         S(std::forward<StorageOrContainer>(s)));
 }
 
 /// histogram factory from vector-like with default storage
-template <typename Iterable, typename = detail::requires_axis_vector<Iterable>>
+template <typename Iterable, typename = detail::requires_sequence_of_any_axis<Iterable>>
 auto make_histogram(Iterable&& c) {
   return make_histogram_with(default_storage(), std::forward<Iterable>(c));
 }
 
 /// histogram factory from vector-like with default storage
-template <typename Iterable, typename = detail::requires_axis_vector<Iterable>>
+template <typename Iterable, typename = detail::requires_sequence_of_any_axis<Iterable>>
 auto make_weighted_histogram(Iterable&& c) {
   return make_histogram_with(weight_storage(), std::forward<Iterable>(c));
 }
@@ -70,7 +73,7 @@ auto make_weighted_histogram(Iterable&& c) {
 template <typename StorageOrContainer, typename Iterator,
           typename = detail::requires_iterator<Iterator>>
 auto make_histogram_with(StorageOrContainer&& s, Iterator begin, Iterator end) {
-  using T = detail::iterator_value_type<Iterator>;
+  using T = detail::unqual<decltype(*begin)>;
   return make_histogram_with(std::forward<StorageOrContainer>(s),
                              boost::container::vector<T>(begin, end));
 }
diff --git a/include/boost/histogram/make_profile.hpp b/include/boost/histogram/make_profile.hpp
index 3c09e013..5a2b199d 100644
--- a/include/boost/histogram/make_profile.hpp
+++ b/include/boost/histogram/make_profile.hpp
@@ -9,7 +9,7 @@
 
 #include <boost/histogram/accumulators/mean.hpp>
 #include <boost/histogram/accumulators/weighted_mean.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/histogram/make_histogram.hpp>
 
 namespace boost {
@@ -30,13 +30,13 @@ auto make_weighted_profile(T&& axis, Ts&&... axes) {
 }
 
 /// profile factory from vector-like
-template <typename Iterable, typename = detail::requires_axis_vector<Iterable>>
+template <typename Iterable, typename = detail::requires_sequence_of_any_axis<Iterable>>
 auto make_profile(Iterable&& c) {
   return make_histogram_with(profile_storage(), std::forward<Iterable>(c));
 }
 
 /// profile factory from vector-like with weighted mean
-template <typename Iterable, typename = detail::requires_axis_vector<Iterable>>
+template <typename Iterable, typename = detail::requires_sequence_of_any_axis<Iterable>>
 auto make_weighted_profile(Iterable&& c) {
   return make_histogram_with(weighted_profile_storage(), std::forward<Iterable>(c));
 }
diff --git a/include/boost/histogram/ostream_operators.hpp b/include/boost/histogram/ostream_operators.hpp
index c88b6cda..7c9ad1a1 100644
--- a/include/boost/histogram/ostream_operators.hpp
+++ b/include/boost/histogram/ostream_operators.hpp
@@ -7,18 +7,32 @@
 #ifndef BOOST_HISTOGRAM_OSTREAM_OPERATORS_HPP
 #define BOOST_HISTOGRAM_OSTREAM_OPERATORS_HPP
 
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <ostream>
 
 namespace boost {
 namespace histogram {
 
+namespace detail {
+template <class OStream, class T>
+void grid_ostream_impl(OStream& os, const char* prefix, const T& t) {
+  os << prefix << "(";
+  t.for_each_axis([&](const auto& a) { os << "\n  " << a << ","; });
+  os << (t.rank() ? "\n)" : ")");
+}
+} // namespace detail
+
+template <typename CharT, typename Traits, typename A, typename S>
+std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os,
+                                              const grid<A, S>& h) {
+  detail::grid_ostream_impl(os, "grid", h);
+  return os;
+}
+
 template <typename CharT, typename Traits, typename A, typename S>
 std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os,
                                               const histogram<A, S>& h) {
-  os << "histogram(";
-  h.for_each_axis([&](const auto& a) { os << "\n  " << a << ","; });
-  os << (h.rank() ? "\n)" : ")");
+  detail::grid_ostream_impl(os, "histogram", h);
   return os;
 }
 
diff --git a/include/boost/histogram/storage_adaptor.hpp b/include/boost/histogram/storage_adaptor.hpp
index 424104a0..6a3ee04e 100644
--- a/include/boost/histogram/storage_adaptor.hpp
+++ b/include/boost/histogram/storage_adaptor.hpp
@@ -11,7 +11,7 @@
 #include <boost/assert.hpp>
 #include <boost/histogram/detail/cat.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 #include <boost/histogram/weight.hpp>
 #include <boost/iterator/iterator_facade.hpp>
 #include <boost/throw_exception.hpp>
diff --git a/test/algorithm_project_test.cpp b/test/algorithm_project_test.cpp
index 6636bed5..edde9c12 100644
--- a/test/algorithm_project_test.cpp
+++ b/test/algorithm_project_test.cpp
@@ -8,7 +8,6 @@
 #include <boost/histogram/algorithm/project.hpp>
 #include <boost/histogram/algorithm/sum.hpp>
 #include <boost/histogram/axis/integer.hpp>
-#include <boost/histogram/axis/ostream_operators.hpp>
 #include <boost/histogram/literals.hpp>
 #include <vector>
 #include "utility_histogram.hpp"
diff --git a/test/algorithm_reduce_test.cpp b/test/algorithm_reduce_test.cpp
index 99ab0e2c..cc85cc52 100644
--- a/test/algorithm_reduce_test.cpp
+++ b/test/algorithm_reduce_test.cpp
@@ -8,11 +8,9 @@
 #include <boost/histogram/algorithm/reduce.hpp>
 #include <boost/histogram/algorithm/sum.hpp>
 #include <boost/histogram/axis/integer.hpp>
-#include <boost/histogram/axis/ostream_operators.hpp>
 #include <boost/histogram/axis/regular.hpp>
 #include <boost/histogram/axis/variable.hpp>
 #include <boost/histogram/indexed.hpp>
-#include <boost/histogram/ostream_operators.hpp>
 #include <boost/histogram/unsafe_access.hpp>
 #include <vector>
 #include "utility_histogram.hpp"
diff --git a/test/detail_test.cpp b/test/detail_test.cpp
index 1887b704..8085ca7d 100644
--- a/test/detail_test.cpp
+++ b/test/detail_test.cpp
@@ -13,7 +13,6 @@
 #include <boost/histogram/axis/variable.hpp>
 #include <boost/histogram/detail/axes.hpp>
 #include <boost/histogram/detail/cat.hpp>
-#include <boost/histogram/detail/is_set.hpp>
 #include <boost/histogram/literals.hpp>
 #include <tuple>
 #include <vector>
@@ -35,43 +34,6 @@ int main() {
     BOOST_TEST_EQ(decltype(213_c)::value, 213);
   }
 
-  // sequence equality
-  {
-    enum { A, B, C };
-    std::vector<axis::variant<axis::regular<>, axis::variable<>, axis::category<>,
-                              axis::integer<>>>
-        std_vector1 = {axis::regular<>{2, -1, 1}, axis::variable<>{-1, 0, 1},
-                       axis::category<>{A, B, C}};
-
-    std::vector<axis::variant<axis::regular<>, axis::variable<>, axis::category<>>>
-        std_vector2 = {axis::regular<>{2, -1, 1}, axis::variable<>{-1, 0, 1},
-                       axis::category<>{{A, B, C}}};
-
-    std::vector<axis::variant<axis::regular<>, axis::variable<>>> std_vector3 = {
-        axis::variable<>{-1, 0, 1}, axis::regular<>{2, -1, 1}};
-
-    std::vector<axis::variant<axis::variable<>, axis::regular<>>> std_vector4 = {
-        axis::regular<>{2, -1, 1}, axis::variable<>{-1, 0, 1}};
-
-    BOOST_TEST(detail::axes_equal(std_vector1, std_vector2));
-    BOOST_TEST_NOT(detail::axes_equal(std_vector2, std_vector3));
-    BOOST_TEST_NOT(detail::axes_equal(std_vector3, std_vector4));
-
-    auto tuple1 = std::make_tuple(axis::regular<>{2, -1, 1}, axis::variable<>{-1, 0, 1},
-                                  axis::category<>{{A, B, C}});
-
-    auto tuple2 = std::make_tuple(axis::regular<>{2, -1, 1}, axis::variable<>{-1, 0, 1},
-                                  axis::category<>{{A, B}});
-
-    auto tuple3 = std::make_tuple(axis::regular<>{2, -1, 1}, axis::variable<>{-1, 0, 1});
-
-    BOOST_TEST(detail::axes_equal(std_vector1, tuple1));
-    BOOST_TEST(detail::axes_equal(tuple1, std_vector1));
-    BOOST_TEST_NOT(detail::axes_equal(tuple1, tuple2));
-    BOOST_TEST_NOT(detail::axes_equal(tuple2, tuple3));
-    BOOST_TEST_NOT(detail::axes_equal(std_vector3, tuple3));
-  }
-
   // axes_size
   {
     std::tuple<int, int> a;
@@ -82,60 +44,102 @@ int main() {
     BOOST_TEST_EQ(detail::axes_size(c), 4);
   }
 
+  // sequence equality
+  {
+    using R = axis::regular<>;
+    using I = axis::integer<>;
+    using V = axis::variable<>;
+    auto r = R(2, -1, 1);
+    auto i = I(-1, 1);
+    auto v = V{-1, 0, 1};
+
+    std::vector<axis::variant<R, I, V>> v1 = {r, i};
+    std::vector<axis::variant<R, I>> v2 = {r, i};
+    std::vector<axis::variant<R, I>> v3 = {i, r};
+    std::vector<axis::variant<I, R>> v4 = {r, i};
+    std::vector<axis::variant<R, I>> v5 = {r, r};
+    std::vector<R> v6 = {r, r};
+
+    BOOST_TEST(detail::axes_equal(v1, v2));
+    BOOST_TEST(detail::axes_equal(v1, v4));
+    BOOST_TEST(detail::axes_equal(v5, v6));
+    BOOST_TEST_NOT(detail::axes_equal(v1, v3));
+    BOOST_TEST_NOT(detail::axes_equal(v2, v3));
+    BOOST_TEST_NOT(detail::axes_equal(v3, v4));
+    BOOST_TEST_NOT(detail::axes_equal(v1, v5));
+
+    auto t1 = std::make_tuple(r, i);
+    auto t2 = std::make_tuple(i, r);
+    auto t3 = std::make_tuple(v, i);
+    auto t4 = std::make_tuple(r, r);
+
+    BOOST_TEST(detail::axes_equal(t1, v1));
+    BOOST_TEST(detail::axes_equal(t1, v2));
+    BOOST_TEST(detail::axes_equal(t1, v4));
+    BOOST_TEST(detail::axes_equal(v1, t1));
+    BOOST_TEST(detail::axes_equal(v2, t1));
+    BOOST_TEST(detail::axes_equal(v4, t1));
+    BOOST_TEST(detail::axes_equal(t2, v3));
+    BOOST_TEST(detail::axes_equal(v3, t2));
+    BOOST_TEST(detail::axes_equal(t4, v5));
+    BOOST_TEST(detail::axes_equal(t4, v6));
+    BOOST_TEST_NOT(detail::axes_equal(t1, t2));
+    BOOST_TEST_NOT(detail::axes_equal(t2, t3));
+    BOOST_TEST_NOT(detail::axes_equal(t1, v3));
+    BOOST_TEST_NOT(detail::axes_equal(t1, v3));
+    BOOST_TEST_NOT(detail::axes_equal(t3, v1));
+    BOOST_TEST_NOT(detail::axes_equal(t3, v2));
+    BOOST_TEST_NOT(detail::axes_equal(t3, v3));
+    BOOST_TEST_NOT(detail::axes_equal(t3, v4));
+  }
+
   // sequence assign
   {
-    enum { A, B, C, D };
-    std::vector<axis::variant<axis::regular<>, axis::variable<>, axis::category<>,
-                              axis::integer<>>>
-        std_vector1 = {axis::regular<>{2, -1, 1}, axis::variable<>{-1, 0, 1},
-                       axis::category<>{A, B, C}};
+    using R = axis::regular<>;
+    using I = axis::integer<>;
+    using V = axis::variable<>;
+    auto r = R(2, -1, 1);
+    auto i = I(-1, 1);
+    auto v = V{-1, 0, 1};
 
-    std::vector<axis::variant<axis::regular<>, axis::variable<>, axis::category<>>>
-        std_vector2 = {axis::regular<>{2, -2, 2}, axis::variable<>{-2, 0, 2},
-                       axis::category<>{A, B}};
+    std::vector<axis::variant<R, V, I>> v1 = {r, i};
+    std::vector<axis::variant<I, R>> v2;
+    std::vector<R> v3 = {r, r};
 
-    detail::axes_assign(std_vector2, std_vector1);
-    BOOST_TEST(detail::axes_equal(std_vector2, std_vector1));
+    BOOST_TEST_NOT(detail::axes_equal(v2, v1));
+    detail::axes_assign(v2, v1);
+    BOOST_TEST(detail::axes_equal(v2, v1));
+    detail::axes_assign(v2, v3);
+    BOOST_TEST(detail::axes_equal(v2, v3));
 
-    auto tuple1 = std::make_tuple(axis::regular<>{2, -3, 3}, axis::variable<>{-3, 0, 3},
-                                  axis::category<>{A, B, C, D});
+    auto t1 = std::make_tuple(r);
+    detail::axes_assign(v3, t1);
+    BOOST_TEST(detail::axes_equal(v3, t1));
 
-    detail::axes_assign(tuple1, std_vector1);
-    BOOST_TEST(detail::axes_equal(tuple1, std_vector1));
+    auto t2 = std::make_tuple(r, i);
+    detail::axes_assign(v2, t2);
+    BOOST_TEST(detail::axes_equal(v2, t2));
 
-    decltype(std_vector1) std_vector3;
-    BOOST_TEST_NOT(detail::axes_equal(std_vector3, tuple1));
-    detail::axes_assign(std_vector3, tuple1);
-    BOOST_TEST(detail::axes_equal(std_vector3, tuple1));
-
-    auto tuple2 = std::make_tuple(axis::regular<>{2, -1, 1}, axis::variable<>{-1, 0, 1},
-                                  axis::category<>{A, B});
-
-    detail::axes_assign(tuple2, tuple1);
-    BOOST_TEST(detail::axes_equal(tuple2, tuple1));
+    auto t3 = std::make_tuple(R{3, -1, 1}, i);
+    BOOST_TEST_NOT(detail::axes_equal(t2, t3));
+    detail::axes_assign(t2, t3);
+    BOOST_TEST(detail::axes_equal(t2, t3));
   }
 
   // make_sub_axes
   {
+    using I = axis::integer<>;
     using boost::mp11::mp_list;
-    axis::integer<> a0(0, 1), a1(1, 2), a2(2, 3);
-    auto axes = std::make_tuple(a0, a1, a2);
-    BOOST_TEST_EQ(detail::make_sub_axes(axes, 0_c), std::make_tuple(a0));
-    BOOST_TEST_EQ(detail::make_sub_axes(axes, 1_c), std::make_tuple(a1));
-    BOOST_TEST_EQ(detail::make_sub_axes(axes, 2_c), std::make_tuple(a2));
-    BOOST_TEST_EQ(detail::make_sub_axes(axes, 0_c, 1_c), std::make_tuple(a0, a1));
-    BOOST_TEST_EQ(detail::make_sub_axes(axes, 0_c, 2_c), std::make_tuple(a0, a2));
-    BOOST_TEST_EQ(detail::make_sub_axes(axes, 1_c, 2_c), std::make_tuple(a1, a2));
-    BOOST_TEST_EQ(detail::make_sub_axes(axes, 0_c, 1_c, 2_c),
-                  std::make_tuple(a0, a1, a2));
-  }
 
-  // is_set
-  {
-    std::vector<int> v({3, 1, 2});
-    BOOST_TEST(detail::is_set(v.begin(), v.end()));
-    v = {3, 1, 2, 1};
-    BOOST_TEST_NOT(detail::is_set(v.begin(), v.end()));
+    I a0(0, 1), a1(1, 2), a2(2, 3);
+    auto ax = std::make_tuple(a0, a1, a2);
+    BOOST_TEST_EQ(detail::make_sub_axes(ax, 0_c), std::make_tuple(a0));
+    BOOST_TEST_EQ(detail::make_sub_axes(ax, 1_c), std::make_tuple(a1));
+    BOOST_TEST_EQ(detail::make_sub_axes(ax, 2_c), std::make_tuple(a2));
+    BOOST_TEST_EQ(detail::make_sub_axes(ax, 0_c, 1_c), std::make_tuple(a0, a1));
+    BOOST_TEST_EQ(detail::make_sub_axes(ax, 0_c, 2_c), std::make_tuple(a0, a2));
+    BOOST_TEST_EQ(detail::make_sub_axes(ax, 1_c, 2_c), std::make_tuple(a1, a2));
+    BOOST_TEST_EQ(detail::make_sub_axes(ax, 0_c, 1_c, 2_c), std::make_tuple(a0, a1, a2));
   }
 
   return boost::report_errors();
diff --git a/test/histogram_dynamic_test.cpp b/test/histogram_dynamic_test.cpp
index 8a1491f7..ed279e37 100644
--- a/test/histogram_dynamic_test.cpp
+++ b/test/histogram_dynamic_test.cpp
@@ -9,12 +9,10 @@
 #include <boost/histogram/adaptive_storage.hpp>
 #include <boost/histogram/axis/category.hpp>
 #include <boost/histogram/axis/integer.hpp>
-#include <boost/histogram/axis/ostream_operators.hpp>
 #include <boost/histogram/axis/regular.hpp>
 #include <boost/histogram/axis/variant.hpp>
 #include <boost/histogram/literals.hpp>
 #include <boost/histogram/make_histogram.hpp>
-#include <boost/histogram/ostream_operators.hpp>
 #include <boost/histogram/storage_adaptor.hpp>
 #include <cstdlib>
 #include <limits>
@@ -31,7 +29,7 @@ using namespace boost::histogram::literals; // to get _c suffix
 int main() {
   // special stuff that only works with dynamic_tag
 
-  // init
+  // init with vector of axis and vector of axis::variant
   {
     auto v = std::vector<axis::variant<axis::regular<>, axis::integer<>>>();
     v.emplace_back(axis::regular<>(4, -1, 1));
diff --git a/test/histogram_mixed_test.cpp b/test/histogram_mixed_test.cpp
index 08a5259b..2cdf263b 100644
--- a/test/histogram_mixed_test.cpp
+++ b/test/histogram_mixed_test.cpp
@@ -7,11 +7,9 @@
 #include <boost/core/lightweight_test.hpp>
 #include <boost/histogram/adaptive_storage.hpp>
 #include <boost/histogram/axis/integer.hpp>
-#include <boost/histogram/axis/ostream_operators.hpp>
 #include <boost/histogram/axis/regular.hpp>
 #include <boost/histogram/histogram.hpp>
 #include <boost/histogram/literals.hpp>
-#include <boost/histogram/ostream_operators.hpp>
 #include <boost/histogram/storage_adaptor.hpp>
 #include "utility_histogram.hpp"
 
@@ -48,9 +46,9 @@ void run_tests() {
 
   // copy_assign
   {
-    auto a = make(T1{}, axis::regular<>{3, 0, 3}, axis::integer<>(0, 2));
+    auto a = make(T1{}, axis::regular<>{3, 0, 3}, axis::integer<>{0, 2});
     auto b = make_s(T2{}, std::vector<unsigned>(), axis::regular<>{3, 0, 3},
-                    axis::integer<>(0, 2));
+                    axis::integer<>{0, 2});
     a(1, 1);
     BOOST_TEST_NE(a, b);
     b = a;
diff --git a/test/histogram_operators_test.cpp b/test/histogram_operators_test.cpp
new file mode 100644
index 00000000..16e158fd
--- /dev/null
+++ b/test/histogram_operators_test.cpp
@@ -0,0 +1,57 @@
+// Copyright 2015-2018 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 <boost/core/lightweight_test.hpp>
+#include <boost/core/lightweight_test_trait.hpp>
+#include <boost/histogram.hpp>
+#include <type_traits>
+#include "utility_histogram.hpp"
+
+using namespace boost::histogram;
+
+template <typename Tag>
+void run_tests() {
+  auto a = make(Tag(), axis::integer<>(0, 2));
+  auto b = a;
+  a(0);
+  b(1);
+  auto c = a + b;
+  BOOST_TEST_EQ(c.at(0), 1);
+  BOOST_TEST_EQ(c.at(1), 1);
+  c += b;
+  BOOST_TEST_EQ(c.at(0), 1);
+  BOOST_TEST_EQ(c.at(1), 2);
+  auto d = a + b + c;
+  BOOST_TEST_EQ(d.at(0), 2);
+  BOOST_TEST_EQ(d.at(1), 3);
+  BOOST_TEST_TRAIT_TRUE((std::is_same<decltype(d), decltype(a)>));
+  auto e = 3 * a;
+  auto f = b * 2;
+  BOOST_TEST_TRAIT_FALSE((std::is_same<decltype(e), decltype(a)>));
+  BOOST_TEST_TRAIT_FALSE((std::is_same<decltype(f), decltype(a)>));
+  BOOST_TEST_EQ(e.at(0), 3);
+  BOOST_TEST_EQ(e.at(1), 0);
+  BOOST_TEST_EQ(f.at(0), 0);
+  BOOST_TEST_EQ(f.at(1), 2);
+  auto r = 1.0 * a;
+  r += b;
+  r += e;
+  BOOST_TEST_EQ(r.at(0), 4);
+  BOOST_TEST_EQ(r.at(1), 1);
+  BOOST_TEST_EQ(r, a + b + 3 * a);
+  auto s = r / 4;
+  r /= 4;
+  BOOST_TEST_EQ(r.at(0), 1);
+  BOOST_TEST_EQ(r.at(1), 0.25);
+  BOOST_TEST_EQ(r, s);
+}
+
+int main() {
+  run_tests<static_tag>();
+  run_tests<dynamic_tag>();
+
+  return boost::report_errors();
+}
diff --git a/test/histogram_serialization_test.cpp b/test/histogram_serialization_test.cpp
index 6cf85085..47d5d1c4 100644
--- a/test/histogram_serialization_test.cpp
+++ b/test/histogram_serialization_test.cpp
@@ -8,8 +8,6 @@
 #include <boost/archive/text_oarchive.hpp>
 #include <boost/core/lightweight_test.hpp>
 #include <boost/histogram.hpp>
-#include <boost/histogram/axis/ostream_operators.hpp>
-#include <boost/histogram/ostream_operators.hpp>
 #include <boost/histogram/serialization.hpp>
 #include <sstream>
 #include <string>
diff --git a/test/histogram_test.cpp b/test/histogram_test.cpp
index eb0c5231..b0483619 100644
--- a/test/histogram_test.cpp
+++ b/test/histogram_test.cpp
@@ -5,20 +5,7 @@
 // or copy at http://www.boost.org/LICENSE_1_0.txt)
 
 #include <boost/core/lightweight_test.hpp>
-#include <boost/histogram/accumulators/mean.hpp>
-#include <boost/histogram/accumulators/ostream_operators.hpp>
-#include <boost/histogram/accumulators/weighted_mean.hpp>
-#include <boost/histogram/accumulators/weighted_sum.hpp>
-#include <boost/histogram/adaptive_storage.hpp>
-#include <boost/histogram/algorithm/sum.hpp>
-#include <boost/histogram/axis.hpp>
-#include <boost/histogram/axis/ostream_operators.hpp>
-#include <boost/histogram/histogram.hpp>
-#include <boost/histogram/literals.hpp>
-#include <boost/histogram/ostream_operators.hpp>
-#include <boost/histogram/sample.hpp>
-#include <boost/histogram/storage_adaptor.hpp>
-#include <boost/histogram/weight.hpp>
+#include <boost/histogram.hpp>
 #include <sstream>
 #include <stdexcept>
 #include <tuple>
@@ -539,40 +526,6 @@ void run_tests() {
     BOOST_TEST_EQ(a[4], 3);
   }
 
-  // operators
-  {
-    auto a = make(Tag(), axis::integer<>(0, 3));
-    auto b = a;
-    a(0);
-    b(1);
-    auto c = a + b;
-    BOOST_TEST_EQ(c.at(0), 1);
-    BOOST_TEST_EQ(c.at(1), 1);
-    c += b;
-    BOOST_TEST_EQ(c.at(0), 1);
-    BOOST_TEST_EQ(c.at(1), 2);
-    auto d = a + b + c;
-    BOOST_TEST_EQ(d.at(0), 2);
-    BOOST_TEST_EQ(d.at(1), 3);
-    auto e = 3 * a;
-    auto f = b * 2;
-    BOOST_TEST_EQ(e.at(0), 3);
-    BOOST_TEST_EQ(e.at(1), 0);
-    BOOST_TEST_EQ(f.at(0), 0);
-    BOOST_TEST_EQ(f.at(1), 2);
-    auto r = a;
-    r += b;
-    r += e;
-    BOOST_TEST_EQ(r.at(0), 4);
-    BOOST_TEST_EQ(r.at(1), 1);
-    BOOST_TEST_EQ(r, a + b + 3 * a);
-    auto s = r / 4;
-    r /= 4;
-    BOOST_TEST_EQ(r.at(0), 1);
-    BOOST_TEST_EQ(r.at(1), 0.25);
-    BOOST_TEST_EQ(r, s);
-  }
-
   // histogram_ostream
   {
     auto a = make(Tag(), axis::regular<>(3, -1, 1, "r"), axis::integer<>(0, 2, "i"));
diff --git a/test/indexed_test.cpp b/test/indexed_test.cpp
index 835d8aaa..36ba85c3 100644
--- a/test/indexed_test.cpp
+++ b/test/indexed_test.cpp
@@ -6,7 +6,6 @@
 
 #include <boost/core/lightweight_test.hpp>
 #include <boost/histogram/axis/integer.hpp>
-#include <boost/histogram/axis/ostream_operators.hpp>
 #include <boost/histogram/axis/variable.hpp>
 #include <boost/histogram/histogram.hpp>
 #include <boost/histogram/indexed.hpp>
diff --git a/test/meta_test.cpp b/test/meta_test.cpp
index c972a375..3616fbf5 100644
--- a/test/meta_test.cpp
+++ b/test/meta_test.cpp
@@ -323,24 +323,6 @@ int main() {
     BOOST_TEST_TRAIT_TRUE((std::is_same<mp_last<L>, long>));
   }
 
-  // container_value_type
-  {
-    using T1 = std::vector<int>;
-    using U1 = container_value_type<T1>;
-    using T2 = const std::vector<const int>&;
-    using U2 = container_value_type<T2>;
-    BOOST_TEST_TRAIT_TRUE((std::is_same<U1, int>));
-    BOOST_TEST_TRAIT_TRUE((std::is_same<U2, const int>));
-  }
-
-  // iterator_value_type
-  {
-    using T1 = const char*;
-    using T2 = std::iterator<std::random_access_iterator_tag, int>;
-    BOOST_TEST_TRAIT_TRUE((std::is_same<iterator_value_type<T1>, char>));
-    BOOST_TEST_TRAIT_TRUE((std::is_same<iterator_value_type<T2>, int>));
-  }
-
   // args_type
   {
     struct Foo {
@@ -388,21 +370,20 @@ int main() {
     BOOST_TEST_EQ(b, std::make_tuple(2, 3));
   }
 
-  // is_axis_vector
+  // is_sequence_of_axis
   {
     using A = std::vector<bh::axis::regular<>>;
     using B = std::vector<bh::axis::variant<bh::axis::regular<>>>;
     using C = std::vector<int>;
-    using D = const std::vector<bh::axis::variant<bh::axis::integer<>>>;
-    using E = const std::vector<bh::axis::variant<bh::axis::integer<>>>&;
     auto v = std::vector<bh::axis::variant<bh::axis::regular<>, bh::axis::integer<>>>();
-    BOOST_TEST_TRAIT_TRUE((is_axis_vector<A>));
-    BOOST_TEST_TRAIT_TRUE((is_axis_vector<B>));
-    BOOST_TEST_TRAIT_FALSE((is_axis_vector<C>));
-    BOOST_TEST_TRAIT_TRUE((is_axis_vector<D>));
-    BOOST_TEST_TRAIT_TRUE((is_axis_vector<E>));
-    BOOST_TEST_TRAIT_TRUE((is_axis_vector<decltype(v)>));
-    BOOST_TEST_TRAIT_TRUE((is_axis_vector<decltype(std::move(v))>));
+    BOOST_TEST_TRAIT_TRUE((is_sequence_of_any_axis<A>));
+    BOOST_TEST_TRAIT_TRUE((is_sequence_of_axis<A>));
+    BOOST_TEST_TRAIT_FALSE((is_sequence_of_axis_variant<A>));
+    BOOST_TEST_TRAIT_TRUE((is_sequence_of_any_axis<B>));
+    BOOST_TEST_TRAIT_TRUE((is_sequence_of_axis_variant<B>));
+    BOOST_TEST_TRAIT_FALSE((is_sequence_of_axis<B>));
+    BOOST_TEST_TRAIT_FALSE((is_sequence_of_any_axis<C>));
+    BOOST_TEST_TRAIT_TRUE((is_sequence_of_any_axis<decltype(v)>));
   }
 
   // is_weight
diff --git a/test/unsafe_access_test.cpp b/test/unsafe_access_test.cpp
index 893bcd79..26cb1721 100644
--- a/test/unsafe_access_test.cpp
+++ b/test/unsafe_access_test.cpp
@@ -8,7 +8,6 @@
 #include <boost/histogram/axis/regular.hpp>
 #include <boost/histogram/indexed.hpp>
 #include <boost/histogram/make_histogram.hpp>
-#include <boost/histogram/ostream_operators.hpp>
 #include <boost/histogram/unsafe_access.hpp>
 #include "utility_histogram.hpp"
 
diff --git a/test/utility_axis.hpp b/test/utility_axis.hpp
index 1fdf70a4..02f371a1 100644
--- a/test/utility_axis.hpp
+++ b/test/utility_axis.hpp
@@ -8,7 +8,7 @@
 #define BOOST_HISTOGRAM_TEST_UTILITY_AXIS_HPP
 
 #include <boost/core/lightweight_test.hpp>
-#include <boost/histogram/histogram_fwd.hpp>
+#include <boost/histogram/fwd.hpp>
 
 namespace boost {
 namespace histogram {
diff --git a/test/utility_histogram.hpp b/test/utility_histogram.hpp
index 33e3b163..39fcc37a 100644
--- a/test/utility_histogram.hpp
+++ b/test/utility_histogram.hpp
@@ -7,8 +7,10 @@
 #ifndef BOOST_HISTOGRAM_TEST_UTILITY_HISTOGRAM_HPP
 #define BOOST_HISTOGRAM_TEST_UTILITY_HISTOGRAM_HPP
 
+#include <boost/histogram/axis/ostream_operators.hpp>
 #include <boost/histogram/axis/variant.hpp>
 #include <boost/histogram/make_histogram.hpp>
+#include <boost/histogram/ostream_operators.hpp>
 #include <type_traits>
 #include <vector>