diff --git a/include/boost/histogram/algorithm/sum.hpp b/include/boost/histogram/algorithm/sum.hpp
index 5ff437dd..6a1f57d1 100644
--- a/include/boost/histogram/algorithm/sum.hpp
+++ b/include/boost/histogram/algorithm/sum.hpp
@@ -9,6 +9,7 @@
 
 #include <boost/histogram/accumulators/sum.hpp>
 #include <boost/histogram/fwd.hpp>
+#include <boost/histogram/indexed.hpp>
 #include <boost/mp11/utility.hpp>
 #include <numeric>
 #include <type_traits>
@@ -16,24 +17,50 @@
 namespace boost {
 namespace histogram {
 namespace algorithm {
-/** Compute the sum over all histogram cells, including underflow/overflow bins.
 
-  If the value type of the histogram is an integral or floating point type,
-  boost::accumulators::sum<double> is used to compute the sum, else the original value
-  type is used. Compilation fails, if the value type does not support operator+=.
+/** Compute the sum over all histogram cells (underflow/overflow included by default).
 
-  Return type is double if the value type of the histogram is integral or floating point,
-  and the original value type otherwise.
- */
+  The implementation favors accuracy and protection against overflow over speed. If the
+  value type of the histogram is an integral or floating point type,
+  accumulators::sum<double> is used to compute the sum, else the original value type is
+  used. Compilation fails, if the value type does not support operator+=. The return type
+  is double if the value type of the histogram is integral or floating point, and the
+  original value type otherwise.
+
+  If you need a different trade-off, you can write your own loop or use `std::accumulate`:
+  ```
+  // iterate over all bins
+  auto sum_all = std::accumulate(hist.begin(), hist.end(), 0.0);
+
+  // skip underflow/overflow bins
+  double sum = 0;
+  for (auto&& x : indexed(hist))
+    sum += *x; // dereference accessor
+
+  // or:
+  // auto ind = boost::histogram::indexed(hist);
+  // auto sum = std::accumulate(ind.begin(), ind.end(), 0.0);
+  ```
+
+  @returns accumulator type or double
+
+  @param hist Const reference to the histogram.
+  @param cov  Iterate over all or only inner bins (optional, default: all).
+*/
 template <class A, class S>
-auto sum(const histogram<A, S>& h) {
+auto sum(const histogram<A, S>& hist, const coverage cov = coverage::all) {
   using T = typename histogram<A, S>::value_type;
-  using Sum = mp11::mp_if<std::is_arithmetic<T>, accumulators::sum<double>, T>;
-  Sum sum;
-  for (auto&& x : h) sum += x;
+  using sum_type = mp11::mp_if<std::is_arithmetic<T>, accumulators::sum<double>, T>;
+  sum_type sum;
+  if (cov == coverage::all)
+    for (auto&& x : hist) sum += x;
+  else
+    // sum += x also works if sum_type::operator+=(const sum_type&) exists
+    for (auto&& x : indexed(hist)) sum += *x;
   using R = mp11::mp_if<std::is_arithmetic<T>, double, T>;
   return static_cast<R>(sum);
 }
+
 } // namespace algorithm
 } // namespace histogram
 } // namespace boost
diff --git a/test/algorithm_sum_test.cpp b/test/algorithm_sum_test.cpp
index 53f80ac9..4f4b8095 100644
--- a/test/algorithm_sum_test.cpp
+++ b/test/algorithm_sum_test.cpp
@@ -9,9 +9,9 @@
 #include <boost/histogram/accumulators/weighted_sum.hpp>
 #include <boost/histogram/algorithm/sum.hpp>
 #include <boost/histogram/axis/integer.hpp>
-#include "throw_exception.hpp"
 #include <unordered_map>
 #include <vector>
+#include "throw_exception.hpp"
 #include "utility_histogram.hpp"
 
 using namespace boost::histogram;
@@ -19,35 +19,56 @@ using boost::histogram::algorithm::sum;
 
 template <typename Tag>
 void run_tests() {
-  auto ax = axis::integer<>(0, 100);
+  auto ax = axis::integer<>(0, 10);
 
-  auto h1 = make(Tag(), ax);
-  for (unsigned i = 0; i < 100; ++i) h1(i);
-  BOOST_TEST_EQ(sum(h1), 100);
+  {
+    auto h = make(Tag(), ax);
+    std::fill(h.begin(), h.end(), 1);
+    BOOST_TEST_EQ(sum(h), 12);
+    BOOST_TEST_EQ(sum(h, coverage::inner), 10);
+  }
 
-  auto h2 = make_s(Tag(), std::vector<double>(), ax, ax);
-  for (unsigned i = 0; i < 100; ++i)
-    for (unsigned j = 0; j < 100; ++j) h2(i, j);
-  BOOST_TEST_EQ(sum(h2), 10000);
+  {
+    auto h = make_s(Tag(), std::array<int, 12>(), ax);
+    std::fill(h.begin(), h.end(), 1);
+    BOOST_TEST_EQ(sum(h), 12);
+    BOOST_TEST_EQ(sum(h, coverage::inner), 10);
+  }
 
-  auto h3 = make_s(Tag(), std::array<int, 102>(), ax);
-  for (unsigned i = 0; i < 100; ++i) h3(i);
-  BOOST_TEST_EQ(sum(h3), 100);
+  {
+    auto h = make_s(Tag(), std::unordered_map<std::size_t, int>(), ax);
+    std::fill(h.begin(), h.end(), 1);
+    BOOST_TEST_EQ(sum(h), 12);
+    BOOST_TEST_EQ(sum(h, coverage::inner), 10);
+  }
 
-  auto h4 = make_s(Tag(), std::unordered_map<std::size_t, int>(), ax);
-  for (unsigned i = 0; i < 100; ++i) h4(i);
-  BOOST_TEST_EQ(sum(h4), 100);
+  {
+    auto h = make_s(Tag(), std::vector<double>(), ax, ax);
+    std::fill(h.begin(), h.end(), 1);
+    BOOST_TEST_EQ(sum(h), 12 * 12);
+    BOOST_TEST_EQ(sum(h, coverage::inner), 10 * 10);
+  }
 
-  auto h5 =
-      make_s(Tag(), std::vector<accumulators::weighted_sum<>>(), axis::integer<>(0, 1),
-             axis::integer<int, axis::null_type, axis::option::none_t>(2, 4));
-  h5(weight(2), 0, 2);
-  h5(-1, 2);
-  h5(1, 3);
+  {
+    using W = accumulators::weighted_sum<>;
+    auto h = make_s(Tag(), std::vector<W>(), axis::integer<>(0, 2),
+                    axis::integer<int, axis::null_type, axis::option::none_t>(2, 4));
+    W w(0, 2);
+    for (auto&& x : h) {
+      x = w;
+      w = W(w.value() + 1, 2);
+    }
 
-  const auto v = algorithm::sum(h5);
-  BOOST_TEST_EQ(v.value(), 4);
-  BOOST_TEST_EQ(v.variance(), 6);
+    // x-axis has 4 bins, y-axis has 2 = 8 bins total with 4 inner bins
+
+    const auto v1 = algorithm::sum(h);
+    BOOST_TEST_EQ(v1.value(), 0 + 1 + 2 + 3 + 4 + 5 + 6 + 7);
+    BOOST_TEST_EQ(v1.variance(), 8 * 2);
+
+    const auto v2 = algorithm::sum(h, coverage::inner);
+    BOOST_TEST_EQ(v2.value(), 1 + 2 + 5 + 6);
+    BOOST_TEST_EQ(v2.variance(), 4 * 2);
+  }
 }
 
 int main() {