diff --git a/build/CMakeLists.txt b/build/CMakeLists.txt
index 8701e4a0..abe4b49d 100644
--- a/build/CMakeLists.txt
+++ b/build/CMakeLists.txt
@@ -14,6 +14,8 @@ endif()
 option(BUILD_PYTHON "Build python bindings" ON)
 option(BUILD_NUMPY_SUPPORT "Build numpy support" ON)
 option(BUILD_CHECKS "Build auxillary checks" OFF)
+option(TRACE_ALLOCS "Debug: Trace allocations" OFF)
+mark_as_advanced(TRACE_ALLOCS)
 
 if(${CMAKE_BUILD_TYPE})
   STRING(TOLOWER ${CMAKE_BUILD_TYPE} CMAKE_BUILD_TYPE)
@@ -26,6 +28,10 @@ else()
   add_definitions(-Wall -fpermissive)
 endif()
 
+if(TRACE_ALLOCS)
+  add_definitions(-DBOOST_HISTOGRAM_TRACE_ALLOCS)
+endif()
+
 if(BUILD_PYTHON)
   if(DEFINED PYTHON_VERSION)
     find_package(PythonLibs ${PYTHON_VERSION} EXACT REQUIRED)
diff --git a/doc/guide.qbk b/doc/guide.qbk
index ef8281f0..7b7f62e2 100644
--- a/doc/guide.qbk
+++ b/doc/guide.qbk
@@ -236,13 +236,26 @@ int main() {
 
 Some arithmetic operations are supported for histograms. Histograms are...
 
-* comparable
+* equal comparable
 * addable (adding histograms with non-matching axes is an error)
 * multipliable and dividable by a real number
 
 The operations are commutative, except for division. Dividing a number by a histogram is not allowed.
 
-The operations are illustrated in the next listing.
+Two histograms compare equal, if...
+
+* all axes compare equal, which also checks equality of axis labels
+* all values and variances compare equal
+
+Adding histograms is useful, if you want to parallelize the filling of a histogram over several threads or processes. Fill independent copies of the histogram in worker threads, and then add them all up in the main thread.
+
+Scaling is useful to re-weight histograms before adding them, for those who need to work with weights. Scaling a bin count by a factor `x` has a different effect on value and variance of the count. The value is multiplied by `x`, but the variance is multiplied by `x*x`. This follows from the properties of the variance, as explained in [link histogram.rationale.variance the rationale].
+
+[warning Because of behavior of the variance, adding a histogram to itself is not identical to multiplying the original histogram by two, as far as the variance is concerned.]
+
+[note Scaling a histogram internally converts the bin counters from integers to double counters per bin to store the now different numbers for value and variance, as explained in previous sections.]
+
+Here is an example which show-cases the supported operators.
 
 [c++]``
 #include <boost/histogram.hpp>
@@ -283,20 +296,30 @@ int main() {
     std::cout << (h4 != h5) << " " << (h4 == 5 * h5) << std::endl;
 
     // prints: 1 1
+
+    // beware the special effect of multiplying a histogram on its variance
+    auto h = bh::make_static_histogram(bh::axis::regular<>(2, -1, 1));
+    h.fill(-0.5);
+    std::cout << "value    " << (2 * h).value(0)
+              << " " << (h + h).value(0) << "\n"
+              << "variance " << (2 * h).variance(0)
+              << " " << (h + h).variance(0) << std::endl;
+    // equality operator also checks variances, so the statement is false
+    std::cout << (h + h == 2 * h) << std::endl;
+
+    /* prints:
+        value    2 2
+        variance 4 2
+        0
+    */
 }
 ``
 
-Adding histograms is useful, if you want to parallelize the filling of a histogram over several threads or processes. Fill independent copies of the histogram in worker threads, and then add them all up in the main thread.
-
-Scaling is useful to re-weight histograms before adding them, for those who need to work with weights.
-
-[note Scaling by a factor `x` internally converts the bin counters from integers to double counters per bin, as explained in previous sections. The value counter is multiplied by `x`, the variance counter by `x*x`. This follows from the properties of the variance, as explained in [link histogram.rationale.variance the rationale].]
-
 [endsect]
 
 [section Streaming]
 
-Simple ostream operators are shipped with the library, which are internally used by the Python interface bindings. These may be useful for debugging, since they give a text representation of axis and histogram instances. The streaming operators are not included by default (so that users can implement their own). The following example shows the effect of output streaming.
+Simple ostream operators are shipped with the library, which are internally used by the Python interface bindings. These give text representations of axis and histogram configurations, but do not show the histogram content. For users, the builtin streaming operators may be useful for debugging. The streaming operators are not included by the standard header `#include <boost/histogram.hpp>`, to not collide with user implementations. The following example shows the effect of output streaming.
 
 [c++]``
 #include <boost/histogram.hpp>
@@ -340,7 +363,7 @@ int main() {
 
 [section Serialization]
 
-The library supports serialization via [@boost:/libs/serialization/index.html Boost.Serialization]. The serialization machinery is also used in Python to enable pickling of histograms. The streaming code is not included by default, so that the library can be used without having Boost.Serialization installed.
+The library supports serialization via [@boost:/libs/serialization/index.html Boost.Serialization]. The serialization machinery is used in the Python module to enable pickling of histograms. The streaming code is not included by default, so that the library can be used without having Boost.Serialization installed.
 
 [c++]``
 #include <boost/histogram.hpp>
@@ -398,13 +421,13 @@ The C++ histogram has Python-bindings, so you can create and fill histograms in
 
 To access the Python interface of Boost.Histogram, you need to active the compilation of the corresponding module, which you can then import via `import histogram` in Python.
 
-The Python interface generally mimics the C++ interface. All methods on the C++ side have an equivalent on the Python side. Typical C++ idioms, however, have been replaced with equivalent Python idioms.
+The Python interface generally mimics the C++ interface. All methods on the C++ side have an equivalent on the Python side. Typical C++ idioms have been replaced with equivalent Python idioms:
 
 * methods which take no argument and just return a value are represented by properties in Python
 * axis objects support the sequence protocol (e.g. len(x) instead of x.size)
 * optional values are passed via keyword arguments
 
-The documentation of the Python interface can be explored from within the Python interpreter (type `help()`, then `histogram`).
+The documentation of the Python interface is best explored from within the Python interpreter (type `help()`, then `histogram`).
 
 Here is an example of the Python module in action.
 
@@ -547,7 +570,7 @@ h2 = hg.histogram(h1) # creates copy
 h1.fill(-0.5)
 h2.fill(0.5)
 
-# arithmetic operators
+# arithmetic operators (see performance note below)
 h3 = h1 + h2
 h4 = h3 * 2
 
@@ -558,6 +581,10 @@ print h4.value(0), h4.value(1)
 ``
 [endsect]
 
+[note Python has no concept of rvalue references and therefore cannot avoid creating temporary objects in arithmetic operations like C++ can. A statement `h3 = h1 + h2` is equivalent to `tmp = hg.histogram(h1); tmp += h2; h3 = hg.histogram(tmp)`, which amounts to two allocations, one in the first and one in the last copy-assignment, while only one allocation is really necessary.
+
+To avoid creating superfluous temporaries, write your code entirely with operators `+=` and `*=`, this is always possible. In the example above, the optimised code would be: `h3 = hg.histogram(h1); h3 += h2`. It is actually quite puzzling why Python does not do such simple optimisations by itself.]
+
 [section Mixing Python and C++]
 
 It is a efficient workflow to create and configure histograms in Python and then pass them to some C++ code which fills them at maximum speed. You rarely need to change the way the histogram is filled, but you likely want to iterate the range and binning of the axis after seeing the data. With [@boost:/libs/python/index.html Boost.Python] it is easy to set this up.
diff --git a/include/boost/histogram/axis.hpp b/include/boost/histogram/axis.hpp
index 72241569..9c9cdaa0 100644
--- a/include/boost/histogram/axis.hpp
+++ b/include/boost/histogram/axis.hpp
@@ -199,6 +199,8 @@ struct pow {
   bool operator==(const pow &other) const noexcept {
     return value == other.value;
   }
+private:
+  friend ::boost::serialization::access;
   template <class Archive> void serialize(Archive &, unsigned);
 };
 } // namespace transform
@@ -570,7 +572,8 @@ public:
   /// Returns the value for the bin index.
   bin_type operator[](int idx) const {
     auto it = map_->right.find(idx);
-    BOOST_ASSERT_MSG(it != map_->right.end(), "category index out of range");
+    if (it == map_->right.end())
+      throw std::out_of_range("category index out of range");
     return it->second;
   }
 
diff --git a/include/boost/histogram/detail/weight.hpp b/include/boost/histogram/detail/weight_counter.hpp
similarity index 53%
rename from include/boost/histogram/detail/weight.hpp
rename to include/boost/histogram/detail/weight_counter.hpp
index b9cf6e78..50592c9e 100644
--- a/include/boost/histogram/detail/weight.hpp
+++ b/include/boost/histogram/detail/weight_counter.hpp
@@ -11,38 +11,39 @@ namespace boost {
 namespace histogram {
 namespace detail {
 
-/// Used by nstore to hold a sum of weighted counts and a variance estimate
-struct weight {
+/// Double counter which holds a sum of weights and a sum of squared weights
+struct weight_counter {
   double w, w2;
-  weight() = default;
-  weight(const weight &) = default;
-  weight(weight &&) = default;
-  weight &operator=(const weight &) = default;
-  weight &operator=(weight &&) = default;
+  weight_counter() = default;
+  weight_counter(const weight_counter &) = default;
+  weight_counter(weight_counter &&) = default;
+  weight_counter &operator=(const weight_counter &) = default;
+  weight_counter &operator=(weight_counter &&) = default;
 
-  weight(double value, double variance) : w(value), w2(variance) {}
+  weight_counter(double value, double variance) : w(value), w2(variance) {}
 
-  weight &operator+=(const weight &rhs) {
-    w += rhs.w;
-    w2 += rhs.w2;
-    return *this;
-  }
-  weight &operator++() {
+  weight_counter &operator++() {
     ++w;
     ++w2;
     return *this;
   }
 
-  weight &operator*=(const double x) {
-    w *= x;
-    w2 *= x;
+  weight_counter &operator+=(const weight_counter &rhs) {
+    w += rhs.w;
+    w2 += rhs.w2;
     return *this;
   }
 
-  bool operator==(const weight &rhs) const {
+  weight_counter &operator*=(const double x) {
+    w *= x;
+    w2 *= x*x;
+    return *this;
+  }
+
+  bool operator==(const weight_counter &rhs) const {
     return w == rhs.w && w2 == rhs.w2;
   }
-  bool operator!=(const weight &rhs) const { return !operator==(rhs); }
+  bool operator!=(const weight_counter &rhs) const { return !operator==(rhs); }
   template <typename T> bool operator==(const T &rhs) const {
     return w == static_cast<double>(rhs) && w2 == static_cast<double>(rhs);
   }
@@ -50,33 +51,33 @@ struct weight {
     return !operator==(rhs);
   }
 
-  weight &add_weight(double t) {
+  weight_counter &operator+=(double t) {
     w += t;
     w2 += t * t;
     return *this;
   }
 
   template <typename T>
-  explicit weight(const T &t)
+  explicit weight_counter(const T &t)
       : w(static_cast<double>(t)), w2(static_cast<double>(t)) {}
 
-  template <typename T> weight &operator=(const T &t) {
+  template <typename T> weight_counter &operator=(const T &t) {
     w = static_cast<double>(t);
     w2 = static_cast<double>(t);
     return *this;
   }
-  template <typename T> weight &operator+=(const T &t) {
+  template <typename T> weight_counter &operator+=(const T &t) {
     w += static_cast<double>(t);
     w2 += static_cast<double>(t);
     return *this;
   }
 };
 
-template <typename T> bool operator==(const T &t, const weight &w) {
+template <typename T> bool operator==(const T &t, const weight_counter &w) {
   return w == t;
 }
 
-template <typename T> bool operator!=(const T &t, const weight &w) {
+template <typename T> bool operator!=(const T &t, const weight_counter &w) {
   return !(w == t);
 }
 } // namespace detail
diff --git a/include/boost/histogram/histogram_impl_dynamic.hpp b/include/boost/histogram/histogram_impl_dynamic.hpp
index 5fcd06a2..8256e60e 100644
--- a/include/boost/histogram/histogram_impl_dynamic.hpp
+++ b/include/boost/histogram/histogram_impl_dynamic.hpp
@@ -8,7 +8,6 @@
 #define _BOOST_HISTOGRAM_HISTOGRAM_DYNAMIC_IMPL_HPP_
 
 #include <algorithm>
-#include <boost/assert.hpp>
 #include <boost/config.hpp>
 #include <boost/histogram/axis.hpp>
 #include <boost/histogram/detail/axis_visitor.hpp>
@@ -133,17 +132,16 @@ public:
     static_assert(
         (n_count::value + n_weight::value) <= 1,
         "arguments may contain at most one instance of type count or weight");
-    BOOST_ASSERT_MSG(sizeof...(args) ==
-                         (dim() + n_count::value + n_weight::value),
-                     "number of arguments does not match histogram dimension");
+    if (dim() != sizeof...(args) - n_count::value - n_weight::value)
+      throw std::invalid_argument("fill arguments does not match histogram dimension");
     fill_impl(mpl::int_<(n_count::value + 2 * n_weight::value)>(),
               std::forward<Args>(args)...);
   }
 
   template <typename Iterator, typename = detail::is_iterator<Iterator>>
   void fill(Iterator begin, Iterator end) {
-    BOOST_ASSERT_MSG(std::distance(begin, end) == dim(),
-                     "number of arguments does not match histogram dimension");
+    if(dim() != std::distance(begin, end))
+      throw std::invalid_argument("fill iterator range does not match histogram dimension");
     std::size_t idx = 0, stride = 1;
     xlin_iter(idx, stride, begin);
     if (stride) {
@@ -153,8 +151,8 @@ public:
 
   template <typename Iterator, typename = detail::is_iterator<Iterator>>
   void fill(Iterator begin, Iterator end, const count n) {
-    BOOST_ASSERT_MSG(std::distance(begin, end) == dim(),
-                     "number of arguments does not match histogram dimension");
+    if(dim() != std::distance(begin, end))
+      throw std::invalid_argument("fill iterator range does not match histogram dimension");
     std::size_t idx = 0, stride = 1;
     xlin_iter(idx, stride, begin);
     if (stride) {
@@ -164,8 +162,8 @@ public:
 
   template <typename Iterator, typename = detail::is_iterator<Iterator>>
   void fill(Iterator begin, Iterator end, const weight w) {
-    BOOST_ASSERT_MSG(std::distance(begin, end) == dim(),
-                     "number of arguments does not match histogram dimension");
+    if(dim() != std::distance(begin, end))
+      throw std::invalid_argument("fill iterator range does not match histogram dimension");
     std::size_t idx = 0, stride = 1;
     xlin_iter(idx, stride, begin);
     if (stride) {
@@ -174,32 +172,30 @@ public:
   }
 
   template <typename... Indices> value_type value(Indices &&... indices) const {
-    BOOST_ASSERT_MSG(sizeof...(indices) == dim(),
-                     "number of arguments does not match histogram dimension");
+    if(dim() != sizeof...(indices))
+      throw std::invalid_argument("value arguments does not match histogram dimension");
     std::size_t idx = 0, stride = 1;
     lin<0>(idx, stride, std::forward<Indices>(indices)...);
-    if (stride == 0) {
+    if (stride == 0)
       throw std::out_of_range("invalid index");
-    }
     return storage_.value(idx);
   }
 
   template <typename Iterator, typename = detail::is_iterator<Iterator>>
   value_type value(Iterator begin, Iterator end) const {
-    BOOST_ASSERT_MSG(std::distance(begin, end) == dim(),
-                     "number of arguments does not match histogram dimension");
+    if(dim() != std::distance(begin, end))
+      throw std::invalid_argument("value iterator range does not match histogram dimension");
     std::size_t idx = 0, stride = 1;
     lin_iter(idx, stride, begin);
-    if (stride == 0) {
+    if (stride == 0)
       throw std::out_of_range("invalid index");
-    }
     return storage_.value(idx);
   }
 
   template <typename... Indices>
   value_type variance(Indices &&... indices) const {
-    BOOST_ASSERT_MSG(sizeof...(indices) == dim(),
-                     "number of arguments does not match histogram dimension");
+    if(dim() != sizeof...(indices))
+      throw std::invalid_argument("variance arguments does not match histogram dimension");
     std::size_t idx = 0, stride = 1;
     lin<0>(idx, stride, std::forward<Indices>(indices)...);
     if (stride == 0) {
@@ -210,8 +206,8 @@ public:
 
   template <typename Iterator, typename = detail::is_iterator<Iterator>>
   value_type variance(Iterator begin, Iterator end) const {
-    BOOST_ASSERT_MSG(std::distance(begin, end) == dim(),
-                     "number of arguments does not match histogram dimension");
+    if(dim() != std::distance(begin, end))
+      throw std::invalid_argument("variance iterator range does not match histogram dimension");
     std::size_t idx = 0, stride = 1;
     lin_iter(idx, stride, begin);
     if (stride == 0) {
@@ -241,13 +237,15 @@ public:
 
   /// Return axis \a i
   any_axis_type &axis(unsigned i = 0) {
-    BOOST_ASSERT_MSG(i < dim(), "axis index out of range");
+    if(i >= dim())
+      throw std::out_of_range("axis index out of range");
     return axes_[i];
   }
 
   /// Return axis \a i (const version)
   const any_axis_type &axis(unsigned i = 0) const {
-    BOOST_ASSERT_MSG(i < dim(), "axis index out of range");
+    if(i >= dim())
+      throw std::out_of_range("axis index out of range");
     return axes_[i];
   }
 
diff --git a/include/boost/histogram/histogram_impl_static.hpp b/include/boost/histogram/histogram_impl_static.hpp
index 6e6908e1..211ea2ae 100644
--- a/include/boost/histogram/histogram_impl_static.hpp
+++ b/include/boost/histogram/histogram_impl_static.hpp
@@ -152,9 +152,8 @@ public:
   /// Sum of all counts in the histogram
   double sum() const noexcept {
     double result = 0.0;
-    for (std::size_t i = 0, n = storage_.size(); i < n; ++i) {
+    for (std::size_t i = 0, n = storage_.size(); i < n; ++i)
       result += storage_.value(i);
-    }
     return result;
   }
 
@@ -210,9 +209,8 @@ private:
   inline void fill_impl(mpl::int_<0>, Args &&... args) {
     std::size_t idx = 0, stride = 1;
     xlin<0>(idx, stride, std::forward<Args>(args)...);
-    if (stride) {
+    if (stride)
       storage_.increase(idx);
-    }
   }
 
   template <typename... Args>
@@ -220,9 +218,8 @@ private:
     std::size_t idx = 0, stride = 1;
     unsigned n = 0;
     xlin_n<0>(idx, stride, n, std::forward<Args>(args)...);
-    if (stride) {
+    if (stride)
       storage_.add(idx, n);
-    }
   }
 
   template <typename... Args>
@@ -230,9 +227,8 @@ private:
     std::size_t idx = 0, stride = 1;
     double w = 0.0;
     xlin_w<0>(idx, stride, w, std::forward<Args>(args)...);
-    if (stride) {
+    if (stride)
       storage_.weighted_increase(idx, w);
-    }
   }
 
   template <unsigned D> inline void lin(std::size_t &, std::size_t &) const {}
diff --git a/include/boost/histogram/serialization.hpp b/include/boost/histogram/serialization.hpp
index 21f2d1ab..e39b3a3c 100644
--- a/include/boost/histogram/serialization.hpp
+++ b/include/boost/histogram/serialization.hpp
@@ -10,7 +10,7 @@
 #include <boost/fusion/algorithm/iteration/for_each.hpp>
 #include <boost/fusion/include/for_each.hpp>
 #include <boost/histogram/detail/utility.hpp>
-#include <boost/histogram/detail/weight.hpp>
+#include <boost/histogram/detail/weight_counter.hpp>
 #include <boost/histogram/histogram.hpp>
 #include <boost/histogram/storage/adaptive_storage.hpp>
 #include <boost/histogram/storage/array_storage.hpp>
@@ -30,7 +30,7 @@ namespace histogram {
 namespace detail {
 
 template <class Archive>
-void serialize(Archive &ar, weight &wt, unsigned /* version */) {
+void serialize(Archive &ar, weight_counter &wt, unsigned /* version */) {
   ar &wt.w;
   ar &wt.w2;
 }
@@ -52,67 +52,67 @@ void serialize(Archive &ar, array_storage<Container> &store,
 template <class Archive>
 void adaptive_storage::serialize(Archive &ar, unsigned /* version */) {
   using detail::array;
-  std::size_t size = this->size();
+  auto size = this->size();
   ar &size;
   if (Archive::is_loading::value) {
-    unsigned tid = 0;
+    auto tid = 0u;
     ar &tid;
-    if (tid == 0) {
+    if (tid == 0u) {
       buffer_ = detail::array<void>(size);
-    } else if (tid == 1) {
+    } else if (tid == 1u) {
       array<uint8_t> a(size);
       ar &serialization::make_array(a.begin(), size);
       buffer_ = std::move(a);
-    } else if (tid == 2) {
+    } else if (tid == 2u) {
       array<uint16_t> a(size);
       ar &serialization::make_array(a.begin(), size);
       buffer_ = std::move(a);
-    } else if (tid == 3) {
+    } else if (tid == 3u) {
       array<uint32_t> a(size);
       ar &serialization::make_array(a.begin(), size);
       buffer_ = std::move(a);
-    } else if (tid == 4) {
+    } else if (tid == 4u) {
       array<uint64_t> a(size);
       ar &serialization::make_array(a.begin(), size);
       buffer_ = std::move(a);
-    } else if (tid == 5) {
+    } else if (tid == 5u) {
       array<detail::mp_int> a(size);
       ar &serialization::make_array(a.begin(), size);
       buffer_ = std::move(a);
-    } else if (tid == 6) {
-      array<detail::weight> a(size);
+    } else if (tid == 6u) {
+      array<detail::weight_counter> a(size);
       ar &serialization::make_array(a.begin(), size);
       buffer_ = std::move(a);
     }
   } else {
-    unsigned tid = 0;
+    auto tid = 0u;
     if (get<array<void>>(&buffer_)) {
-      tid = 0;
+      tid = 0u;
       ar &tid;
-    } else if (array<uint8_t> *a = get<array<uint8_t>>(&buffer_)) {
-      tid = 1;
+    } else if (auto *a = get<array<uint8_t>>(&buffer_)) {
+      tid = 1u;
       ar &tid;
       ar &serialization::make_array(a->begin(), size);
-    } else if (array<uint16_t> *a = get<array<uint16_t>>(&buffer_)) {
-      tid = 2;
+    } else if (auto *a = get<array<uint16_t>>(&buffer_)) {
+      tid = 2u;
       ar &tid;
       ar &serialization::make_array(a->begin(), size);
-    } else if (array<uint32_t> *a = get<array<uint32_t>>(&buffer_)) {
-      tid = 3;
+    } else if (auto *a = get<array<uint32_t>>(&buffer_)) {
+      tid = 3u;
       ar &tid;
       ar &serialization::make_array(a->begin(), size);
-    } else if (array<uint64_t> *a = get<array<uint64_t>>(&buffer_)) {
-      tid = 4;
+    } else if (auto *a = get<array<uint64_t>>(&buffer_)) {
+      tid = 4u;
       ar &tid;
       ar &serialization::make_array(a->begin(), size);
-    } else if (array<detail::mp_int> *a =
+    } else if (auto *a =
                    get<array<detail::mp_int>>(&buffer_)) {
-      tid = 5;
+      tid = 5u;
       ar &tid;
       ar &serialization::make_array(a->begin(), size);
-    } else if (array<detail::weight> *a =
-                   get<array<detail::weight>>(&buffer_)) {
-      tid = 6;
+    } else if (auto *a =
+                   get<array<detail::weight_counter>>(&buffer_)) {
+      tid = 6u;
       ar &tid;
       ar &serialization::make_array(a->begin(), size);
     }
diff --git a/include/boost/histogram/storage/adaptive_storage.hpp b/include/boost/histogram/storage/adaptive_storage.hpp
index 946b31f7..c48f9479 100644
--- a/include/boost/histogram/storage/adaptive_storage.hpp
+++ b/include/boost/histogram/storage/adaptive_storage.hpp
@@ -10,13 +10,18 @@
 #include <algorithm>
 #include <boost/cstdint.hpp>
 #include <boost/histogram/detail/meta.hpp>
-#include <boost/histogram/detail/weight.hpp>
+#include <boost/histogram/detail/weight_counter.hpp>
 #include <boost/mpl/int.hpp>
 #include <boost/multiprecision/cpp_int.hpp>
 #include <boost/variant.hpp>
 #include <limits>
 #include <memory>
 #include <type_traits>
+#ifdef BOOST_HISTOGRAM_TRACE_ALLOCS
+#include <iostream>
+#include <boost/core/typeinfo.hpp>
+#endif
+
 
 // forward declaration for serialization
 namespace boost {
@@ -39,6 +44,15 @@ namespace detail {
 
 using mp_int = multiprecision::cpp_int;
 
+template <typename T>
+inline T* alloc(std::size_t s) {
+#ifdef BOOST_HISTOGRAM_TRACE_ALLOCS
+  boost::core::typeinfo const & ti = BOOST_CORE_TYPEID(T);
+  std::cerr << "alloc " << boost::core::demangled_name( ti ) << "[" << s << "]" << std::endl;
+#endif
+  return new T[s];
+}
+
 class array_base {
 public:
   explicit array_base(const std::size_t s) : size(s) {}
@@ -58,18 +72,18 @@ public:
 
 template <typename T> class array : public array_base {
 public:
-  explicit array(const std::size_t s) : array_base(s), ptr(new T[s]) {
+  explicit array(const std::size_t s) : array_base(s), ptr(alloc<T>(s)) {
     std::fill(begin(), end(), T(0));
   }
   array() = default;
-  array(const array &rhs) : array_base(rhs), ptr(new T[rhs.size]) {
+  array(const array &rhs) : array_base(rhs), ptr(alloc<T>(rhs.size)) {
     std::copy(rhs.begin(), rhs.end(), begin());
   }
   array &operator=(const array &rhs) {
     if (this != &rhs) {
       if (size != rhs.size) {
         size = rhs.size;
-        ptr.reset(new T[size]);
+        ptr.reset(alloc<T>(size));
       }
       std::copy(rhs.begin(), rhs.end(), begin());
     }
@@ -91,7 +105,7 @@ public:
   template <typename U>
   array(const array<U> &rhs,
         std::size_t nmax = std::numeric_limits<std::size_t>::max())
-      : array_base(rhs), ptr(new T[rhs.size]) {
+      : array_base(rhs), ptr(alloc<T>(rhs.size)) {
     std::copy(rhs.begin(), rhs.begin() + std::min(nmax, size), begin());
   }
 
@@ -114,7 +128,7 @@ public:
 
 using any_array =
     variant<array<void>, array<uint8_t>, array<uint16_t>, array<uint32_t>,
-            array<uint64_t>, array<mp_int>, array<weight>>;
+            array<uint64_t>, array<mp_int>, array<weight_counter>>;
 
 template <typename T> struct next_type;
 template <> struct next_type<uint8_t> { using type = uint16_t; };
@@ -180,7 +194,7 @@ template <typename RHS> struct assign_visitor : public static_visitor<void> {
 
   void operator()(array<mp_int> &lhs) const { lhs[idx].assign(rhs); }
 
-  void operator()(array<weight> &lhs) const { lhs[idx] = rhs; }
+  void operator()(array<weight_counter> &lhs) const { lhs[idx] = rhs; }
 };
 
 struct increase_visitor : public static_visitor<void> {
@@ -204,7 +218,7 @@ struct increase_visitor : public static_visitor<void> {
 
   void operator()(array<mp_int> &lhs) const { ++lhs[idx]; }
 
-  void operator()(array<weight> &lhs) const { ++lhs[idx]; }
+  void operator()(array<weight_counter> &lhs) const { ++lhs[idx]; }
 };
 
 struct wincrease_visitor : public static_visitor<void> {
@@ -215,18 +229,18 @@ struct wincrease_visitor : public static_visitor<void> {
       : lhs_any(l), idx(i), rhs(r) {}
 
   template <typename T> void operator()(array<T> &lhs) const {
-    array<weight> a(lhs);
-    a[idx].add_weight(rhs);
+    array<weight_counter> a(lhs);
+    a[idx] += rhs;
     lhs_any = std::move(a);
   }
 
   void operator()(array<void> &lhs) const {
-    array<weight> a(lhs.size);
-    a[idx].add_weight(rhs);
+    array<weight_counter> a(lhs.size);
+    a[idx] += rhs;
     lhs_any = std::move(a);
   }
 
-  void operator()(array<weight> &lhs) const { lhs[idx].add_weight(rhs); }
+  void operator()(array<weight_counter> &lhs) const { lhs[idx] += rhs; }
 };
 
 struct value_visitor : public static_visitor<double> {
@@ -239,7 +253,7 @@ struct value_visitor : public static_visitor<double> {
 
   double operator()(const array<void> & /*b*/) const { return 0; }
 
-  double operator()(const array<weight> &b) const { return b[idx].w; }
+  double operator()(const array<weight_counter> &b) const { return b[idx].w; }
 };
 
 struct variance_visitor : public static_visitor<double> {
@@ -252,7 +266,7 @@ struct variance_visitor : public static_visitor<double> {
 
   double operator()(const array<void> & /*b*/) const { return 0; }
 
-  double operator()(const array<weight> &b) const { return b[idx].w2; }
+  double operator()(const array<weight_counter> &b) const { return b[idx].w2; }
 };
 
 template <typename RHS> struct radd_visitor : public static_visitor<void> {
@@ -280,27 +294,27 @@ template <typename RHS> struct radd_visitor : public static_visitor<void> {
     lhs[idx] += static_cast<mp_int>(rhs);
   }
 
-  void operator()(array<weight> &lhs) const { lhs[idx] += rhs; }
+  void operator()(array<weight_counter> &lhs) const { lhs[idx] += rhs; }
 };
 
-template <> struct radd_visitor<weight> : public static_visitor<void> {
+template <> struct radd_visitor<weight_counter> : public static_visitor<void> {
   any_array &lhs_any;
   const std::size_t idx;
-  const weight &rhs;
-  radd_visitor(any_array &l, const std::size_t i, const weight &r)
+  const weight_counter &rhs;
+  radd_visitor(any_array &l, const std::size_t i, const weight_counter &r)
       : lhs_any(l), idx(i), rhs(r) {}
 
   template <typename T> void operator()(array<T> &lhs) const {
-    lhs_any = array<weight>(lhs);
-    operator()(get<array<weight>>(lhs_any));
+    lhs_any = array<weight_counter>(lhs);
+    operator()(get<array<weight_counter>>(lhs_any));
   }
 
   void operator()(array<void> &lhs) const {
-    lhs_any = array<weight>(lhs.size);
-    operator()(get<array<weight>>(lhs_any));
+    lhs_any = array<weight_counter>(lhs.size);
+    operator()(get<array<weight_counter>>(lhs_any));
   }
 
-  void operator()(array<weight> &lhs) const { lhs[idx] += rhs; }
+  void operator()(array<weight_counter> &lhs) const { lhs[idx] += rhs; }
 };
 
 // precondition: both arrays must have same size and may not be identical
@@ -319,11 +333,11 @@ struct rmul_visitor : public static_visitor<void> {
   const double x;
   rmul_visitor(any_array &l, const double v) : lhs_any(l), x(v) {}
   template <typename T> void operator()(array<T> &lhs) const {
-    lhs_any = array<weight>(lhs);
-    operator()(get<array<weight>>(lhs_any));
+    lhs_any = array<weight_counter>(lhs);
+    operator()(get<array<weight_counter>>(lhs_any));
   }
   void operator()(array<void> &) const {}
-  void operator()(array<weight> &lhs) const {
+  void operator()(array<weight_counter> &lhs) const {
     for (auto i = 0ul; i != lhs.size; ++i)
       lhs[i] *= x;
   }
@@ -412,14 +426,14 @@ public:
     if (value == variance) {
       apply_visitor(detail::radd_visitor<T>(buffer_, i, value), buffer_);
     } else {
-      apply_visitor(detail::radd_visitor<detail::weight>(
-                        buffer_, i, detail::weight(value, variance)),
+      apply_visitor(detail::radd_visitor<detail::weight_counter>(
+                        buffer_, i, detail::weight_counter(value, variance)),
                     buffer_);
     }
   }
 
-  void weighted_increase(std::size_t i, const double weight) {
-    apply_visitor(detail::wincrease_visitor(buffer_, i, weight), buffer_);
+  void weighted_increase(std::size_t i, const double weight_counter) {
+    apply_visitor(detail::wincrease_visitor(buffer_, i, weight_counter), buffer_);
   }
 
   value_type value(std::size_t i) const {
diff --git a/src/python/histogram.cpp b/src/python/histogram.cpp
index 4cbb4467..e5bef8da 100644
--- a/src/python/histogram.cpp
+++ b/src/python/histogram.cpp
@@ -38,7 +38,7 @@ namespace python {
 class access {
 public:
   using mp_int = histogram::detail::mp_int;
-  using weight = histogram::detail::weight;
+  using weight_counter = histogram::detail::weight_counter;
   template <typename T>
   using array = histogram::detail::array<T>;
 
@@ -58,7 +58,7 @@ public:
       strides.append(sizeof(double));
       return dtype_typestr<double>();
     }
-    str operator()(const array<weight>& /*unused*/) const {
+    str operator()(const array<weight_counter>& /*unused*/) const {
       strides.append(sizeof(double));
       strides.append(strides[-1] * 2);
       shapes.append(2);
diff --git a/test/adaptive_storage_test.cpp b/test/adaptive_storage_test.cpp
index 4dc8a862..498172d0 100644
--- a/test/adaptive_storage_test.cpp
+++ b/test/adaptive_storage_test.cpp
@@ -32,7 +32,7 @@ template <> adaptive_storage prepare<void>(unsigned n) {
   return s;
 }
 
-template <> adaptive_storage prepare<detail::weight>(unsigned n) {
+template <> adaptive_storage prepare<detail::weight_counter>(unsigned n) {
   adaptive_storage s(n);
   s.weighted_increase(0, 1.0);
   return s;
@@ -309,7 +309,7 @@ int main() {
 
   // copy
   {
-    copy_impl<detail::weight>();
+    copy_impl<detail::weight_counter>();
     copy_impl<void>();
     copy_impl<uint8_t>();
     copy_impl<uint16_t>();
@@ -326,7 +326,7 @@ int main() {
     equal_impl<uint32_t>();
     equal_impl<uint64_t>();
     equal_impl<detail::mp_int>();
-    equal_impl<detail::weight>();
+    equal_impl<detail::weight_counter>();
   }
 
   // increase_and_grow
@@ -379,20 +379,20 @@ int main() {
     adaptive_storage a(2);
     a.increase(0);
     a *= 3;
-    BOOST_TEST_EQ(a.value(0), 3.0);
-    BOOST_TEST_EQ(a.variance(0), 3.0);
-    BOOST_TEST_EQ(a.value(1), 0.0);
-    BOOST_TEST_EQ(a.variance(1), 0.0);
-    a.add(1, 2.0, 5.0);
-    BOOST_TEST_EQ(a.value(0), 3.0);
-    BOOST_TEST_EQ(a.variance(0), 3.0);
-    BOOST_TEST_EQ(a.value(1), 2.0);
-    BOOST_TEST_EQ(a.variance(1), 5.0);
+    BOOST_TEST_EQ(a.value(0), 3);
+    BOOST_TEST_EQ(a.variance(0), 9);
+    BOOST_TEST_EQ(a.value(1), 0);
+    BOOST_TEST_EQ(a.variance(1), 0);
+    a.add(1, 2, 5);
+    BOOST_TEST_EQ(a.value(0), 3);
+    BOOST_TEST_EQ(a.variance(0), 9);
+    BOOST_TEST_EQ(a.value(1), 2);
+    BOOST_TEST_EQ(a.variance(1), 5);
     a *= 3;
-    BOOST_TEST_EQ(a.value(0), 9.0);
-    BOOST_TEST_EQ(a.variance(0), 9.0);
-    BOOST_TEST_EQ(a.value(1), 6.0);
-    BOOST_TEST_EQ(a.variance(1), 15.0);
+    BOOST_TEST_EQ(a.value(0), 9);
+    BOOST_TEST_EQ(a.variance(0), 81);
+    BOOST_TEST_EQ(a.value(1), 6);
+    BOOST_TEST_EQ(a.variance(1), 45);
   }
 
   // convert_array_storage
@@ -403,7 +403,7 @@ int main() {
     convert_array_storage_impl<uint32_t>();
     convert_array_storage_impl<uint64_t>();
     convert_array_storage_impl<detail::mp_int>();
-    convert_array_storage_impl<detail::weight>();
+    convert_array_storage_impl<detail::weight_counter>();
   }
 
   // serialization_test
@@ -414,7 +414,7 @@ int main() {
     serialization_impl<uint32_t>();
     serialization_impl<uint64_t>();
     serialization_impl<detail::mp_int>();
-    serialization_impl<detail::weight>();
+    serialization_impl<detail::weight_counter>();
   }
 
   return boost::report_errors();
diff --git a/test/detail_test.cpp b/test/detail_test.cpp
index 309e5284..e4c9dbd0 100644
--- a/test/detail_test.cpp
+++ b/test/detail_test.cpp
@@ -7,23 +7,41 @@
 #include <boost/core/lightweight_test.hpp>
 #include <boost/histogram/detail/axis_visitor.hpp>
 #include <boost/histogram/detail/utility.hpp>
-#include <boost/histogram/detail/weight.hpp>
+#include <boost/histogram/detail/weight_counter.hpp>
 #include <boost/variant.hpp>
 #include <cstring>
 #include <sstream>
 using namespace boost::histogram::detail;
 
+namespace boost { namespace histogram { namespace detail {
+std::ostream& operator<<(std::ostream& os, const weight_counter& w) {
+  os << "[ " << w.w << ", " << w.w2 << "]";
+  return os;
+}
+}}}
+
 int main() {
-  // weight
+  // weight_counter
   {
-    BOOST_TEST(weight(0) == weight());
-    weight w(1);
-    BOOST_TEST(w == weight(1));
-    BOOST_TEST(w != weight());
-    BOOST_TEST(1 == w);
-    BOOST_TEST(w == 1);
-    BOOST_TEST(2 != w);
-    BOOST_TEST(w != 2);
+    BOOST_TEST_EQ(weight_counter(0), weight_counter());
+    weight_counter w(1);
+    BOOST_TEST_EQ(w, weight_counter(1));
+    BOOST_TEST_NE(w, weight_counter());
+    BOOST_TEST_EQ(1, w);
+    BOOST_TEST_EQ(w, 1);
+    BOOST_TEST_NE(2, w);
+    BOOST_TEST_NE(w, 2);
+    w += 2.0;
+    BOOST_TEST_EQ(w, weight_counter(3, 5));
+    // consistency: a weighted counter increased by weight 1 multiplied
+    // by 2 must be the same as a weighted counter increased by weight 2
+    weight_counter u(0);
+    u += 1.0;
+    u *= 2;
+    BOOST_TEST_EQ(u, weight_counter(2, 4));
+    weight_counter v(0);
+    v += 2.0;
+    BOOST_TEST_EQ(u, v);
   }
 
   // escape0
diff --git a/test/python_suite_test.py b/test/python_suite_test.py
index d80bbbbf..48a5407d 100644
--- a/test/python_suite_test.py
+++ b/test/python_suite_test.py
@@ -591,12 +591,14 @@ class test_histogram(unittest.TestCase):
         self.assertEqual(h.value(0), 2)
         self.assertEqual(h.variance(0), 2)
         self.assertEqual(h.value(1), 0)
-        h *= 2.5
-        self.assertEqual(h.value(0), 5)
-        self.assertEqual(h.variance(0), 5)
+        h *= 2
+        self.assertEqual(h.value(0), 4)
+        self.assertEqual(h.variance(0), 8)
         self.assertEqual(h.value(1), 0)
-        self.assertEqual(h + h, 2 * h)
-        self.assertEqual(h + h, h * 2)
+        self.assertEqual((h + h).value(0), (2 * h).value(0))
+        self.assertEqual((h + h).value(0), (h * 2).value(0))
+        self.assertNotEqual((h + h).variance(0), (2 * h).variance(0))
+        self.assertNotEqual((h + h).variance(0), (h * 2).variance(0))
         h2 = histogram(regular(2, 0, 2))
         with self.assertRaises(RuntimeError):
             h + h2