boost/histogram/unlimited_storage.hpp

   1 // Copyright 2015-2019 Hans Dembinski
   2 // Copyright 2019 Glen Joseph Fernandes (glenjofe@gmail.com)
   3 //
   4 // Distributed under the Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt
   6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #ifndef BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
   9 #define BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
  10
  11 #include <algorithm>
  12 #include <boost/assert.hpp>
  13 #include <boost/config.hpp>
  14 #include <boost/core/alloc_construct.hpp>
  15 #include <boost/core/exchange.hpp>
  16 #include <boost/core/nvp.hpp>
  17 #include <boost/histogram/detail/array_wrapper.hpp>
  18 #include <boost/histogram/detail/iterator_adaptor.hpp>
  19 #include <boost/histogram/detail/large_int.hpp>
  20 #include <boost/histogram/detail/operators.hpp>
  21 #include <boost/histogram/detail/safe_comparison.hpp>
  22 #include <boost/histogram/fwd.hpp>
  23 #include <boost/mp11/algorithm.hpp>
  24 #include <boost/mp11/list.hpp>
  25 #include <boost/mp11/utility.hpp>
  26 #include <cmath>
  27 #include <cstdint>
  28 #include <functional>
  29 #include <iterator>
  30 #include <memory>
  31 #include <type_traits>
  32
  33 namespace boost {
  34 namespace histogram {
  35 namespace detail {
  36
  37 template <class T>
  38 struct is_large_int : std::false_type {};
  39
  40 template <class A>
  41 struct is_large_int<large_int<A>> : std::true_type {};
  42
  43 template <class T, class ReturnType>
  44 using if_arithmetic_or_large_int =
  45     std::enable_if_t<(std::is_arithmetic<T>::value || is_large_int<T>::value),
  46                      ReturnType>;
  47
  48 template <class L, class T>
  49 using next_type = mp11::mp_at_c<L, (mp11::mp_find<L, T>::value + 1)>;
  50
  51 template <class Allocator>
  52 class construct_guard {
  53 public:
  54   using pointer = typename std::allocator_traits<Allocator>::pointer;
  55
  56   construct_guard(Allocator& a, pointer p, std::size_t n) noexcept
  57       : a_(a), p_(p), n_(n) {}
  58
  59   ~construct_guard() {
  60     if (p_) { a_.deallocate(p_, n_); }
  61   }
  62
  63   void release() { p_ = pointer(); }
  64
  65   construct_guard(const construct_guard&) = delete;
  66   construct_guard& operator=(const construct_guard&) = delete;
  67
  68 private:
  69   Allocator& a_;
  70   pointer p_;
  71   std::size_t n_;
  72 };
  73
  74 template <class Allocator>
  75 void* buffer_create(Allocator& a, std::size_t n) {
  76   auto ptr = a.allocate(n); // may throw
  77   static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
  78                 "ptr must be trivially copyable");
  79   construct_guard<Allocator> guard(a, ptr, n);
  80   boost::alloc_construct_n(a, ptr, n);
  81   guard.release();
  82   return static_cast<void*>(ptr);
  83 }
  84
  85 template <class Allocator, class Iterator>
  86 auto buffer_create(Allocator& a, std::size_t n, Iterator iter) {
  87   BOOST_ASSERT(n > 0u);
  88   auto ptr = a.allocate(n); // may throw
  89   static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
  90                 "ptr must be trivially copyable");
  91   construct_guard<Allocator> guard(a, ptr, n);
  92   using T = typename std::allocator_traits<Allocator>::value_type;
  93   struct casting_iterator {
  94     void operator++() noexcept { ++iter_; }
  95     T operator*() noexcept {
  96       return static_cast<T>(*iter_);
  97     } // silence conversion warnings
  98     Iterator iter_;
  99   };
 100   boost::alloc_construct_n(a, ptr, n, casting_iterator{iter});
 101   guard.release();
 102   return ptr;
 103 }
 104
 105 template <class Allocator>
 106 void buffer_destroy(Allocator& a, typename std::allocator_traits<Allocator>::pointer p,
 107                     std::size_t n) {
 108   BOOST_ASSERT(p);
 109   BOOST_ASSERT(n > 0u);
 110   boost::alloc_destroy_n(a, p, n);
 111   a.deallocate(p, n);
 112 }
 113
 114 } // namespace detail
 115
 116 /**
 117   Memory-efficient storage for integral counters which cannot overflow.
 118
 119   This storage provides a no-overflow-guarantee if the counters are incremented with
 120   integer weights. It maintains a contiguous array of elemental counters, one for each
 121   cell. If an operation is requested which would overflow a counter, the array is
 122   replaced with another of a wider integral type, then the operation is executed. The
 123   storage uses integers of 8, 16, 32, 64 bits, and then switches to a multiprecision
 124   integral type, similar to those in
 125   [Boost.Multiprecision](https://www.boost.org/doc/libs/develop/libs/multiprecision/doc/html/index.html).
 126
 127   A scaling operation or adding a floating point number triggers a conversion of the
 128   elemental counters into doubles, which voids the no-overflow-guarantee.
 129 */
 130 template <class Allocator>
 131 class unlimited_storage {
 132   static_assert(
 133       std::is_same<typename std::allocator_traits<Allocator>::pointer,
 134                    typename std::allocator_traits<Allocator>::value_type*>::value,
 135       "unlimited_storage requires allocator with trivial pointer type");
 136   using U8 = std::uint8_t;
 137   using U16 = std::uint16_t;
 138   using U32 = std::uint32_t;
 139   using U64 = std::uint64_t;
 140
 141 public:
 142   static constexpr bool has_threading_support = false;
 143
 144   using allocator_type = Allocator;
 145   using value_type = double;
 146   using large_int = detail::large_int<
 147       typename std::allocator_traits<allocator_type>::template rebind_alloc<U64>>;
 148
 149   struct buffer_type {
 150     // cannot be moved outside of scope of unlimited_storage, large_int is dependent type
 151     using types = mp11::mp_list<U8, U16, U32, U64, large_int, double>;
 152
 153     template <class T>
 154     static constexpr unsigned type_index() noexcept {
 155       return static_cast<unsigned>(mp11::mp_find<types, T>::value);
 156     }
 157
 158     template <class F, class... Ts>
 159     decltype(auto) visit(F&& f, Ts&&... ts) const {
 160       // this is intentionally not a switch, the if-chain is faster in benchmarks
 161       if (type == type_index<U8>())
 162         return f(static_cast<U8*>(ptr), std::forward<Ts>(ts)...);
 163       if (type == type_index<U16>())
 164         return f(static_cast<U16*>(ptr), std::forward<Ts>(ts)...);
 165       if (type == type_index<U32>())
 166         return f(static_cast<U32*>(ptr), std::forward<Ts>(ts)...);
 167       if (type == type_index<U64>())
 168         return f(static_cast<U64*>(ptr), std::forward<Ts>(ts)...);
 169       if (type == type_index<large_int>())
 170         return f(static_cast<large_int*>(ptr), std::forward<Ts>(ts)...);
 171       return f(static_cast<double*>(ptr), std::forward<Ts>(ts)...);
 172     }
 173
 174     buffer_type(const allocator_type& a = {}) : alloc(a) {}
 175
 176     buffer_type(buffer_type&& o) noexcept
 177         : alloc(std::move(o.alloc))
 178         , size(boost::exchange(o.size, 0))
 179         , type(boost::exchange(o.type, 0))
 180         , ptr(boost::exchange(o.ptr, nullptr)) {}
 181
 182     buffer_type& operator=(buffer_type&& o) noexcept {
 183       using std::swap;
 184       swap(alloc, o.alloc);
 185       swap(size, o.size);
 186       swap(type, o.type);
 187       swap(ptr, o.ptr);
 188       return *this;
 189     }
 190
 191     buffer_type(const buffer_type& x) : alloc(x.alloc) {
 192       x.visit([this, n = x.size](const auto* xp) {
 193         using T = std::decay_t<decltype(*xp)>;
 194         this->template make<T>(n, xp);
 195       });
 196     }
 197
 198     buffer_type& operator=(const buffer_type& o) {
 199       *this = buffer_type(o);
 200       return *this;
 201     }
 202
 203     ~buffer_type() noexcept { destroy(); }
 204
 205     void destroy() noexcept {
 206       BOOST_ASSERT((ptr == nullptr) == (size == 0));
 207       if (ptr == nullptr) return;
 208       visit([this](auto* p) {
 209         using T = std::decay_t<decltype(*p)>;
 210         using alloc_type =
 211             typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
 212         alloc_type a(alloc); // rebind allocator
 213         detail::buffer_destroy(a, p, this->size);
 214       });
 215       size = 0;
 216       type = 0;
 217       ptr = nullptr;
 218     }
 219
 220     template <class T>
 221     void make(std::size_t n) {
 222       // note: order of commands is to not leave buffer in invalid state upon throw
 223       destroy();
 224       if (n > 0) {
 225         // rebind allocator
 226         using alloc_type =
 227             typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
 228         alloc_type a(alloc);
 229         ptr = detail::buffer_create(a, n); // may throw
 230       }
 231       size = n;
 232       type = type_index<T>();
 233     }
 234
 235     template <class T, class U>
 236     void make(std::size_t n, U iter) {
 237       // note: iter may be current ptr, so create new buffer before deleting old buffer
 238       void* new_ptr = nullptr;
 239       const auto new_type = type_index<T>();
 240       if (n > 0) {
 241         // rebind allocator
 242         using alloc_type =
 243             typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
 244         alloc_type a(alloc);
 245         new_ptr = detail::buffer_create(a, n, iter); // may throw
 246       }
 247       destroy();
 248       size = n;
 249       type = new_type;
 250       ptr = new_ptr;
 251     }
 252
 253     allocator_type alloc;
 254     std::size_t size = 0;
 255     unsigned type = 0;
 256     mutable void* ptr = nullptr;
 257   };
 258
 259   class reference; // forward declare to make friend of const_reference
 260
 261   /// implementation detail
 262   class const_reference
 263       : detail::partially_ordered<const_reference, const_reference, void> {
 264   public:
 265     const_reference(buffer_type& b, std::size_t i) noexcept : bref_(b), idx_(i) {
 266       BOOST_ASSERT(idx_ < bref_.size);
 267     }
 268
 269     const_reference(const const_reference&) noexcept = default;
 270
 271     // no assignment for const_references
 272     const_reference& operator=(const const_reference&) = delete;
 273     const_reference& operator=(const_reference&&) = delete;
 274
 275     operator double() const noexcept {
 276       return bref_.visit(
 277           [this](const auto* p) { return static_cast<double>(p[this->idx_]); });
 278     }
 279
 280     bool operator<(const const_reference& o) const noexcept {
 281       return apply_binary<detail::safe_less>(o);
 282     }
 283
 284     bool operator==(const const_reference& o) const noexcept {
 285       return apply_binary<detail::safe_equal>(o);
 286     }
 287
 288     template <class U>
 289     detail::if_arithmetic_or_large_int<U, bool> operator<(const U& o) const noexcept {
 290       return apply_binary<detail::safe_less>(o);
 291     }
 292
 293     template <class U>
 294     detail::if_arithmetic_or_large_int<U, bool> operator>(const U& o) const noexcept {
 295       return apply_binary<detail::safe_greater>(o);
 296     }
 297
 298     template <class U>
 299     detail::if_arithmetic_or_large_int<U, bool> operator==(const U& o) const noexcept {
 300       return apply_binary<detail::safe_equal>(o);
 301     }
 302
 303   private:
 304     template <class Binary>
 305     bool apply_binary(const const_reference& x) const noexcept {
 306       return x.bref_.visit([this, ix = x.idx_](const auto* xp) {
 307         return this->apply_binary<Binary>(xp[ix]);
 308       });
 309     }
 310
 311     template <class Binary, class U>
 312     bool apply_binary(const U& x) const noexcept {
 313       return bref_.visit([i = idx_, &x](const auto* p) { return Binary()(p[i], x); });
 314     }
 315
 316   protected:
 317     buffer_type& bref_;
 318     std::size_t idx_;
 319     friend class reference;
 320   };
 321
 322   /// implementation detail
 323   class reference : public const_reference,
 324                     public detail::partially_ordered<reference, reference, void> {
 325   public:
 326     reference(buffer_type& b, std::size_t i) noexcept : const_reference(b, i) {}
 327
 328     // references do copy-construct
 329     reference(const reference& x) noexcept = default;
 330
 331     // references do not rebind, assign through
 332     reference& operator=(const reference& x) {
 333       return operator=(static_cast<const_reference>(x));
 334     }
 335
 336     // references do not rebind, assign through
 337     reference& operator=(const const_reference& x) {
 338       // safe for self-assignment, assigning matching type doesn't invalide buffer
 339       x.bref_.visit([this, ix = x.idx_](const auto* xp) { this->operator=(xp[ix]); });
 340       return *this;
 341     }
 342
 343     template <class U>
 344     detail::if_arithmetic_or_large_int<U, reference&> operator=(const U& x) {
 345       this->bref_.visit([this, &x](auto* p) {
 346         // gcc-8 optimizes the expression `p[this->idx_] = 0` away even at -O0,
 347         // so we merge it into the next line which is properly counted
 348         adder()((p[this->idx_] = 0, p), this->bref_, this->idx_, x);
 349       });
 350       return *this;
 351     }
 352
 353     bool operator<(const reference& o) const noexcept {
 354       return const_reference::operator<(o);
 355     }
 356
 357     bool operator==(const reference& o) const noexcept {
 358       return const_reference::operator==(o);
 359     }
 360
 361     template <class U>
 362     detail::if_arithmetic_or_large_int<U, bool> operator<(const U& o) const noexcept {
 363       return const_reference::operator<(o);
 364     }
 365
 366     template <class U>
 367     detail::if_arithmetic_or_large_int<U, bool> operator>(const U& o) const noexcept {
 368       return const_reference::operator>(o);
 369     }
 370
 371     template <class U>
 372     detail::if_arithmetic_or_large_int<U, bool> operator==(const U& o) const noexcept {
 373       return const_reference::operator==(o);
 374     }
 375
 376     reference& operator+=(const const_reference& x) {
 377       x.bref_.visit([this, ix = x.idx_](const auto* xp) { this->operator+=(xp[ix]); });
 378       return *this;
 379     }
 380
 381     template <class U>
 382     detail::if_arithmetic_or_large_int<U, reference&> operator+=(const U& x) {
 383       this->bref_.visit(adder(), this->bref_, this->idx_, x);
 384       return *this;
 385     }
 386
 387     reference& operator-=(const double x) { return operator+=(-x); }
 388
 389     reference& operator*=(const double x) {
 390       this->bref_.visit(multiplier(), this->bref_, this->idx_, x);
 391       return *this;
 392     }
 393
 394     reference& operator/=(const double x) { return operator*=(1.0 / x); }
 395
 396     reference& operator++() {
 397       this->bref_.visit(incrementor(), this->bref_, this->idx_);
 398       return *this;
 399     }
 400   };
 401
 402 private:
 403   template <class Value, class Reference>
 404   class iterator_impl : public detail::iterator_adaptor<iterator_impl<Value, Reference>,
 405                                                         std::size_t, Reference, Value> {
 406   public:
 407     iterator_impl() = default;
 408     template <class V, class R>
 409     iterator_impl(const iterator_impl<V, R>& it)
 410         : iterator_impl::iterator_adaptor_(it.base()), buffer_(it.buffer_) {}
 411     iterator_impl(buffer_type* b, std::size_t i) noexcept
 412         : iterator_impl::iterator_adaptor_(i), buffer_(b) {}
 413
 414     Reference operator*() const noexcept { return {*buffer_, this->base()}; }
 415
 416     template <class V, class R>
 417     friend class iterator_impl;
 418
 419   private:
 420     mutable buffer_type* buffer_ = nullptr;
 421   };
 422
 423 public:
 424   using const_iterator = iterator_impl<const value_type, const_reference>;
 425   using iterator = iterator_impl<value_type, reference>;
 426
 427   explicit unlimited_storage(const allocator_type& a = {}) : buffer_(a) {}
 428   unlimited_storage(const unlimited_storage&) = default;
 429   unlimited_storage& operator=(const unlimited_storage&) = default;
 430   unlimited_storage(unlimited_storage&&) = default;
 431   unlimited_storage& operator=(unlimited_storage&&) = default;
 432
 433   // TODO
 434   // template <class Allocator>
 435   // unlimited_storage(const unlimited_storage<Allocator>& s)
 436
 437   template <class Iterable, class = detail::requires_iterable<Iterable>>
 438   explicit unlimited_storage(const Iterable& s) {
 439     using std::begin;
 440     using std::end;
 441     auto s_begin = begin(s);
 442     auto s_end = end(s);
 443     using V = typename std::iterator_traits<decltype(begin(s))>::value_type;
 444     constexpr auto ti = buffer_type::template type_index<V>();
 445     constexpr auto nt = mp11::mp_size<typename buffer_type::types>::value;
 446     const std::size_t size = static_cast<std::size_t>(std::distance(s_begin, s_end));
 447 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
 448     if
 449 #else
 450     if constexpr
 451 #endif
 452         (ti < nt)
 453       buffer_.template make<V>(size, s_begin);
 454     else
 455       buffer_.template make<double>(size, s_begin);
 456   }
 457
 458   template <class Iterable, class = detail::requires_iterable<Iterable>>
 459   unlimited_storage& operator=(const Iterable& s) {
 460     *this = unlimited_storage(s);
 461     return *this;
 462   }
 463
 464   allocator_type get_allocator() const { return buffer_.alloc; }
 465
 466   void reset(std::size_t n) { buffer_.template make<U8>(n); }
 467
 468   std::size_t size() const noexcept { return buffer_.size; }
 469
 470   reference operator[](std::size_t i) noexcept { return {buffer_, i}; }
 471   const_reference operator[](std::size_t i) const noexcept { return {buffer_, i}; }
 472
 473   bool operator==(const unlimited_storage& x) const noexcept {
 474     if (size() != x.size()) return false;
 475     return buffer_.visit([&x](const auto* p) {
 476       return x.buffer_.visit([p, n = x.size()](const auto* xp) {
 477         return std::equal(p, p + n, xp, detail::safe_equal{});
 478       });
 479     });
 480   }
 481
 482   template <class Iterable>
 483   bool operator==(const Iterable& iterable) const {
 484     if (size() != iterable.size()) return false;
 485     return buffer_.visit([&iterable](const auto* p) {
 486       return std::equal(p, p + iterable.size(), std::begin(iterable),
 487                         detail::safe_equal{});
 488     });
 489   }
 490
 491   unlimited_storage& operator*=(const double x) {
 492     buffer_.visit(multiplier(), buffer_, x);
 493     return *this;
 494   }
 495
 496   iterator begin() noexcept { return {&buffer_, 0}; }
 497   iterator end() noexcept { return {&buffer_, size()}; }
 498   const_iterator begin() const noexcept { return {&buffer_, 0}; }
 499   const_iterator end() const noexcept { return {&buffer_, size()}; }
 500
 501   /// implementation detail; used by unit tests, not part of generic storage interface
 502   template <class T>
 503   unlimited_storage(std::size_t s, const T* p, const allocator_type& a = {})
 504       : buffer_(std::move(a)) {
 505     buffer_.template make<T>(s, p);
 506   }
 507
 508   template <class Archive>
 509   void serialize(Archive& ar, unsigned /* version */) {
 510     if (Archive::is_loading::value) {
 511       buffer_type tmp(buffer_.alloc);
 512       std::size_t size;
 513       ar& make_nvp("type", tmp.type);
 514       ar& make_nvp("size", size);
 515       tmp.visit([this, size](auto* tp) {
 516         BOOST_ASSERT(tp == nullptr);
 517         using T = std::decay_t<decltype(*tp)>;
 518         buffer_.template make<T>(size);
 519       });
 520     } else {
 521       ar& make_nvp("type", buffer_.type);
 522       ar& make_nvp("size", buffer_.size);
 523     }
 524     buffer_.visit([this, &ar](auto* tp) {
 525       auto w = detail::make_array_wrapper(tp, this->buffer_.size);
 526       ar& make_nvp("buffer", w);
 527     });
 528   }
 529
 530 private:
 531   struct incrementor {
 532     template <class T>
 533     void operator()(T* tp, buffer_type& b, std::size_t i) {
 534       BOOST_ASSERT(tp && i < b.size);
 535       if (!detail::safe_increment(tp[i])) {
 536         using U = detail::next_type<typename buffer_type::types, T>;
 537         b.template make<U>(b.size, tp);
 538         ++static_cast<U*>(b.ptr)[i];
 539       }
 540     }
 541
 542     void operator()(large_int* tp, buffer_type&, std::size_t i) { ++tp[i]; }
 543
 544     void operator()(double* tp, buffer_type&, std::size_t i) { ++tp[i]; }
 545   };
 546
 547   struct adder {
 548     template <class U>
 549     void operator()(double* tp, buffer_type&, std::size_t i, const U& x) {
 550       tp[i] += static_cast<double>(x);
 551     }
 552
 553     void operator()(large_int* tp, buffer_type&, std::size_t i, const large_int& x) {
 554       tp[i] += x; // potentially adding large_int to itself is safe
 555     }
 556
 557     template <class T, class U>
 558     void operator()(T* tp, buffer_type& b, std::size_t i, const U& x) {
 559       is_x_integral(std::is_integral<U>{}, tp, b, i, x);
 560     }
 561
 562     template <class T, class U>
 563     void is_x_integral(std::false_type, T* tp, buffer_type& b, std::size_t i,
 564                        const U& x) {
 565       // x could be reference to buffer we manipulate, make copy before changing buffer
 566       const auto v = static_cast<double>(x);
 567       b.template make<double>(b.size, tp);
 568       operator()(static_cast<double*>(b.ptr), b, i, v);
 569     }
 570
 571     template <class T>
 572     void is_x_integral(std::false_type, T* tp, buffer_type& b, std::size_t i,
 573                        const large_int& x) {
 574       // x could be reference to buffer we manipulate, make copy before changing buffer
 575       const auto v = static_cast<large_int>(x);
 576       b.template make<large_int>(b.size, tp);
 577       operator()(static_cast<large_int*>(b.ptr), b, i, v);
 578     }
 579
 580     template <class T, class U>
 581     void is_x_integral(std::true_type, T* tp, buffer_type& b, std::size_t i, const U& x) {
 582       is_x_unsigned(std::is_unsigned<U>{}, tp, b, i, x);
 583     }
 584
 585     template <class T, class U>
 586     void is_x_unsigned(std::false_type, T* tp, buffer_type& b, std::size_t i,
 587                        const U& x) {
 588       if (x >= 0)
 589         is_x_unsigned(std::true_type{}, tp, b, i, detail::make_unsigned(x));
 590       else
 591         is_x_integral(std::false_type{}, tp, b, i, static_cast<double>(x));
 592     }
 593
 594     template <class T, class U>
 595     void is_x_unsigned(std::true_type, T* tp, buffer_type& b, std::size_t i, const U& x) {
 596       if (detail::safe_radd(tp[i], x)) return;
 597       // x could be reference to buffer we manipulate, need to convert to value
 598       const auto y = x;
 599       using TN = detail::next_type<typename buffer_type::types, T>;
 600       b.template make<TN>(b.size, tp);
 601       is_x_unsigned(std::true_type{}, static_cast<TN*>(b.ptr), b, i, y);
 602     }
 603
 604     template <class U>
 605     void is_x_unsigned(std::true_type, large_int* tp, buffer_type&, std::size_t i,
 606                        const U& x) {
 607       tp[i] += x;
 608     }
 609   };
 610
 611   struct multiplier {
 612     template <class T>
 613     void operator()(T* tp, buffer_type& b, const double x) {
 614       // potential lossy conversion that cannot be avoided
 615       b.template make<double>(b.size, tp);
 616       operator()(static_cast<double*>(b.ptr), b, x);
 617     }
 618
 619     void operator()(double* tp, buffer_type& b, const double x) {
 620       for (auto end = tp + b.size; tp != end; ++tp) *tp *= x;
 621     }
 622
 623     template <class T>
 624     void operator()(T* tp, buffer_type& b, std::size_t i, const double x) {
 625       b.template make<double>(b.size, tp);
 626       operator()(static_cast<double*>(b.ptr), b, i, x);
 627     }
 628
 629     void operator()(double* tp, buffer_type&, std::size_t i, const double x) {
 630       tp[i] *= static_cast<double>(x);
 631     }
 632   };
 633
 634   mutable buffer_type buffer_;
 635   friend struct unsafe_access;
 636 };
 637
 638 } // namespace histogram
 639 } // namespace boost
 640
 641 #endif