Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-12-16 09:53:11

0001 // Copyright 2015-2019 Hans Dembinski
0002 // Copyright 2019 Glen Joseph Fernandes (glenjofe@gmail.com)
0003 //
0004 // Distributed under the Boost Software License, Version 1.0.
0005 // (See accompanying file LICENSE_1_0.txt
0006 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0007 
0008 #ifndef BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
0009 #define BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
0010 
0011 #include <algorithm>
0012 #include <boost/core/alloc_construct.hpp>
0013 #include <boost/core/exchange.hpp>
0014 #include <boost/core/nvp.hpp>
0015 #include <boost/histogram/detail/array_wrapper.hpp>
0016 #include <boost/histogram/detail/iterator_adaptor.hpp>
0017 #include <boost/histogram/detail/large_int.hpp>
0018 #include <boost/histogram/detail/operators.hpp>
0019 #include <boost/histogram/detail/safe_comparison.hpp>
0020 #include <boost/histogram/fwd.hpp>
0021 #include <boost/mp11/algorithm.hpp>
0022 #include <boost/mp11/list.hpp>
0023 #include <boost/mp11/utility.hpp>
0024 #include <cassert>
0025 #include <cmath>
0026 #include <cstdint>
0027 #include <functional>
0028 #include <iterator>
0029 #include <memory>
0030 #include <type_traits>
0031 
0032 namespace boost {
0033 namespace histogram {
0034 namespace detail {
0035 
0036 template <class T>
0037 struct is_large_int : std::false_type {};
0038 
0039 template <class A>
0040 struct is_large_int<large_int<A>> : std::true_type {};
0041 
0042 template <class T, class ReturnType>
0043 using if_arithmetic_or_large_int =
0044     std::enable_if_t<(std::is_arithmetic<T>::value || is_large_int<T>::value),
0045                      ReturnType>;
0046 
0047 template <class L, class T>
0048 using next_type = mp11::mp_at_c<L, (mp11::mp_find<L, T>::value + 1)>;
0049 
0050 template <class Allocator>
0051 class construct_guard {
0052 public:
0053   using pointer = typename std::allocator_traits<Allocator>::pointer;
0054 
0055   construct_guard(Allocator& a, pointer p, std::size_t n) noexcept
0056       : a_(a), p_(p), n_(n) {}
0057 
0058   ~construct_guard() {
0059     if (p_) { a_.deallocate(p_, n_); }
0060   }
0061 
0062   void release() { p_ = pointer(); }
0063 
0064   construct_guard(const construct_guard&) = delete;
0065   construct_guard& operator=(const construct_guard&) = delete;
0066 
0067 private:
0068   Allocator& a_;
0069   pointer p_;
0070   std::size_t n_;
0071 };
0072 
0073 template <class Allocator>
0074 void* buffer_create(Allocator& a, std::size_t n) {
0075   auto ptr = a.allocate(n); // may throw
0076   static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
0077                 "ptr must be trivially copyable");
0078   construct_guard<Allocator> guard(a, ptr, n);
0079   boost::alloc_construct_n(a, ptr, n);
0080   guard.release();
0081   return static_cast<void*>(ptr);
0082 }
0083 
0084 template <class Allocator, class Iterator>
0085 auto buffer_create(Allocator& a, std::size_t n, Iterator iter) {
0086   assert(n > 0u);
0087   auto ptr = a.allocate(n); // may throw
0088   static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
0089                 "ptr must be trivially copyable");
0090   construct_guard<Allocator> guard(a, ptr, n);
0091   using T = typename std::allocator_traits<Allocator>::value_type;
0092   struct casting_iterator {
0093     void operator++() noexcept { ++iter_; }
0094     T operator*() noexcept {
0095       return static_cast<T>(*iter_);
0096     } // silence conversion warnings
0097     Iterator iter_;
0098   };
0099   boost::alloc_construct_n(a, ptr, n, casting_iterator{iter});
0100   guard.release();
0101   return ptr;
0102 }
0103 
0104 template <class Allocator>
0105 void buffer_destroy(Allocator& a, typename std::allocator_traits<Allocator>::pointer p,
0106                     std::size_t n) {
0107   assert(p);
0108   assert(n > 0u);
0109   boost::alloc_destroy_n(a, p, n);
0110   a.deallocate(p, n);
0111 }
0112 
0113 } // namespace detail
0114 
0115 /**
0116   Memory-efficient storage for integral counters which cannot overflow.
0117 
0118   This storage provides a no-overflow-guarantee if the counters are incremented with
0119   integer weights. It maintains a contiguous array of elemental counters, one for each
0120   cell. If an operation is requested which would overflow a counter, the array is
0121   replaced with another of a wider integral type, then the operation is executed. The
0122   storage uses integers of 8, 16, 32, 64 bits, and then switches to a multiprecision
0123   integral type, similar to those in
0124   [Boost.Multiprecision](https://www.boost.org/doc/libs/develop/libs/multiprecision/doc/html/index.html).
0125 
0126   A scaling operation or adding a floating point number triggers a conversion of the
0127   elemental counters into doubles, which voids the no-overflow-guarantee.
0128 */
0129 template <class Allocator>
0130 class unlimited_storage {
0131   static_assert(
0132       std::is_same<typename std::allocator_traits<Allocator>::pointer,
0133                    typename std::allocator_traits<Allocator>::value_type*>::value,
0134       "unlimited_storage requires allocator with trivial pointer type");
0135   using U8 = std::uint8_t;
0136   using U16 = std::uint16_t;
0137   using U32 = std::uint32_t;
0138   using U64 = std::uint64_t;
0139 
0140 public:
0141   static constexpr bool has_threading_support = false;
0142 
0143   using allocator_type = Allocator;
0144   using value_type = double;
0145   using large_int = detail::large_int<
0146       typename std::allocator_traits<allocator_type>::template rebind_alloc<U64>>;
0147 
0148   struct buffer_type {
0149     // cannot be moved outside of scope of unlimited_storage, large_int is dependent type
0150     using types = mp11::mp_list<U8, U16, U32, U64, large_int, double>;
0151 
0152     template <class T>
0153     static constexpr unsigned type_index() noexcept {
0154       return static_cast<unsigned>(mp11::mp_find<types, T>::value);
0155     }
0156 
0157     template <class F, class... Ts>
0158     decltype(auto) visit(F&& f, Ts&&... ts) const {
0159       // this is intentionally not a switch, the if-chain is faster in benchmarks
0160       if (type == type_index<U8>())
0161         return f(static_cast<U8*>(ptr), std::forward<Ts>(ts)...);
0162       if (type == type_index<U16>())
0163         return f(static_cast<U16*>(ptr), std::forward<Ts>(ts)...);
0164       if (type == type_index<U32>())
0165         return f(static_cast<U32*>(ptr), std::forward<Ts>(ts)...);
0166       if (type == type_index<U64>())
0167         return f(static_cast<U64*>(ptr), std::forward<Ts>(ts)...);
0168       if (type == type_index<large_int>())
0169         return f(static_cast<large_int*>(ptr), std::forward<Ts>(ts)...);
0170       return f(static_cast<double*>(ptr), std::forward<Ts>(ts)...);
0171     }
0172 
0173     buffer_type(const allocator_type& a = {}) : alloc(a) {}
0174 
0175     buffer_type(buffer_type&& o) noexcept
0176         : alloc(std::move(o.alloc))
0177         , size(boost::exchange(o.size, 0))
0178         , type(boost::exchange(o.type, 0))
0179         , ptr(boost::exchange(o.ptr, nullptr)) {}
0180 
0181     buffer_type& operator=(buffer_type&& o) noexcept {
0182       using std::swap;
0183       swap(alloc, o.alloc);
0184       swap(size, o.size);
0185       swap(type, o.type);
0186       swap(ptr, o.ptr);
0187       return *this;
0188     }
0189 
0190     buffer_type(const buffer_type& x) : alloc(x.alloc) {
0191       x.visit([this, n = x.size](const auto* xp) {
0192         using T = std::decay_t<decltype(*xp)>;
0193         this->template make<T>(n, xp);
0194       });
0195     }
0196 
0197     buffer_type& operator=(const buffer_type& o) {
0198       *this = buffer_type(o);
0199       return *this;
0200     }
0201 
0202     ~buffer_type() noexcept { destroy(); }
0203 
0204     void destroy() noexcept {
0205       assert((ptr == nullptr) == (size == 0));
0206       if (ptr == nullptr) return;
0207       visit([this](auto* p) {
0208         using T = std::decay_t<decltype(*p)>;
0209         using alloc_type =
0210             typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
0211         alloc_type a(alloc); // rebind allocator
0212         detail::buffer_destroy(a, p, this->size);
0213       });
0214       size = 0;
0215       type = 0;
0216       ptr = nullptr;
0217     }
0218 
0219     template <class T>
0220     void make(std::size_t n) {
0221       // note: order of commands is to not leave buffer in invalid state upon throw
0222       destroy();
0223       if (n > 0) {
0224         // rebind allocator
0225         using alloc_type =
0226             typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
0227         alloc_type a(alloc);
0228         ptr = detail::buffer_create(a, n); // may throw
0229       }
0230       size = n;
0231       type = type_index<T>();
0232     }
0233 
0234     template <class T, class U>
0235     void make(std::size_t n, U iter) {
0236       // note: iter may be current ptr, so create new buffer before deleting old buffer
0237       void* new_ptr = nullptr;
0238       const auto new_type = type_index<T>();
0239       if (n > 0) {
0240         // rebind allocator
0241         using alloc_type =
0242             typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
0243         alloc_type a(alloc);
0244         new_ptr = detail::buffer_create(a, n, iter); // may throw
0245       }
0246       destroy();
0247       size = n;
0248       type = new_type;
0249       ptr = new_ptr;
0250     }
0251 
0252     allocator_type alloc;
0253     std::size_t size = 0;
0254     unsigned type = 0;
0255     mutable void* ptr = nullptr;
0256   };
0257 
0258   class reference; // forward declare to make friend of const_reference
0259 
0260   /// implementation detail
0261   class const_reference
0262       : detail::partially_ordered<const_reference, const_reference, void> {
0263   public:
0264     const_reference(buffer_type& b, std::size_t i) noexcept : bref_(b), idx_(i) {
0265       assert(idx_ < bref_.size);
0266     }
0267 
0268     const_reference(const const_reference&) noexcept = default;
0269 
0270     // no assignment for const_references
0271     const_reference& operator=(const const_reference&) = delete;
0272     const_reference& operator=(const_reference&&) = delete;
0273 
0274     operator double() const noexcept {
0275       return bref_.visit(
0276           [this](const auto* p) { return static_cast<double>(p[this->idx_]); });
0277     }
0278 
0279     bool operator<(const const_reference& o) const noexcept {
0280       return apply_binary<detail::safe_less>(o);
0281     }
0282 
0283     bool operator==(const const_reference& o) const noexcept {
0284       return apply_binary<detail::safe_equal>(o);
0285     }
0286 
0287     template <class U>
0288     detail::if_arithmetic_or_large_int<U, bool> operator<(const U& o) const noexcept {
0289       return apply_binary<detail::safe_less>(o);
0290     }
0291 
0292     template <class U>
0293     detail::if_arithmetic_or_large_int<U, bool> operator>(const U& o) const noexcept {
0294       return apply_binary<detail::safe_greater>(o);
0295     }
0296 
0297     template <class U>
0298     detail::if_arithmetic_or_large_int<U, bool> operator==(const U& o) const noexcept {
0299       return apply_binary<detail::safe_equal>(o);
0300     }
0301 
0302   private:
0303     template <class Binary>
0304     bool apply_binary(const const_reference& x) const noexcept {
0305       return x.bref_.visit([this, ix = x.idx_](const auto* xp) {
0306         return this->apply_binary<Binary>(xp[ix]);
0307       });
0308     }
0309 
0310     template <class Binary, class U>
0311     bool apply_binary(const U& x) const noexcept {
0312       return bref_.visit([i = idx_, &x](const auto* p) { return Binary()(p[i], x); });
0313     }
0314 
0315   protected:
0316     buffer_type& bref_;
0317     std::size_t idx_;
0318     friend class reference;
0319   };
0320 
0321   /// implementation detail
0322   class reference : public const_reference,
0323                     public detail::partially_ordered<reference, reference, void> {
0324   public:
0325     reference(buffer_type& b, std::size_t i) noexcept : const_reference(b, i) {}
0326 
0327     // references do copy-construct
0328     reference(const reference& x) noexcept = default;
0329 
0330     // references do not rebind, assign through
0331     reference& operator=(const reference& x) {
0332       return operator=(static_cast<const_reference>(x));
0333     }
0334 
0335     // references do not rebind, assign through
0336     reference& operator=(const const_reference& x) {
0337       // safe for self-assignment, assigning matching type doesn't invalide buffer
0338       x.bref_.visit([this, ix = x.idx_](const auto* xp) { this->operator=(xp[ix]); });
0339       return *this;
0340     }
0341 
0342     template <class U>
0343     detail::if_arithmetic_or_large_int<U, reference&> operator=(const U& x) {
0344       this->bref_.visit([this, &x](auto* p) {
0345         // gcc-8 optimizes the expression `p[this->idx_] = 0` away even at -O0,
0346         // so we merge it into the next line which is properly counted
0347         adder()((p[this->idx_] = 0, p), this->bref_, this->idx_, x);
0348       });
0349       return *this;
0350     }
0351 
0352     bool operator<(const reference& o) const noexcept {
0353       return const_reference::operator<(o);
0354     }
0355 
0356     bool operator==(const reference& o) const noexcept {
0357       return const_reference::operator==(o);
0358     }
0359 
0360     template <class U>
0361     detail::if_arithmetic_or_large_int<U, bool> operator<(const U& o) const noexcept {
0362       return const_reference::operator<(o);
0363     }
0364 
0365     template <class U>
0366     detail::if_arithmetic_or_large_int<U, bool> operator>(const U& o) const noexcept {
0367       return const_reference::operator>(o);
0368     }
0369 
0370     template <class U>
0371     detail::if_arithmetic_or_large_int<U, bool> operator==(const U& o) const noexcept {
0372       return const_reference::operator==(o);
0373     }
0374 
0375     reference& operator+=(const const_reference& x) {
0376       x.bref_.visit([this, ix = x.idx_](const auto* xp) { this->operator+=(xp[ix]); });
0377       return *this;
0378     }
0379 
0380     template <class U>
0381     detail::if_arithmetic_or_large_int<U, reference&> operator+=(const U& x) {
0382       this->bref_.visit(adder(), this->bref_, this->idx_, x);
0383       return *this;
0384     }
0385 
0386     reference& operator-=(const double x) { return operator+=(-x); }
0387 
0388     reference& operator*=(const double x) {
0389       this->bref_.visit(multiplier(), this->bref_, this->idx_, x);
0390       return *this;
0391     }
0392 
0393     reference& operator/=(const double x) { return operator*=(1.0 / x); }
0394 
0395     reference& operator++() {
0396       this->bref_.visit(incrementor(), this->bref_, this->idx_);
0397       return *this;
0398     }
0399   };
0400 
0401 private:
0402   template <class Value, class Reference>
0403   class iterator_impl : public detail::iterator_adaptor<iterator_impl<Value, Reference>,
0404                                                         std::size_t, Reference, Value> {
0405   public:
0406     iterator_impl() = default;
0407     template <class V, class R>
0408     iterator_impl(const iterator_impl<V, R>& it)
0409         : iterator_impl::iterator_adaptor_(it.base()), buffer_(it.buffer_) {}
0410     iterator_impl(buffer_type* b, std::size_t i) noexcept
0411         : iterator_impl::iterator_adaptor_(i), buffer_(b) {}
0412 
0413     Reference operator*() const noexcept { return {*buffer_, this->base()}; }
0414 
0415     template <class V, class R>
0416     friend class iterator_impl;
0417 
0418   private:
0419     mutable buffer_type* buffer_ = nullptr;
0420   };
0421 
0422 public:
0423   using const_iterator = iterator_impl<const value_type, const_reference>;
0424   using iterator = iterator_impl<value_type, reference>;
0425 
0426   explicit unlimited_storage(const allocator_type& a = {}) : buffer_(a) {}
0427   unlimited_storage(const unlimited_storage&) = default;
0428   unlimited_storage& operator=(const unlimited_storage&) = default;
0429   unlimited_storage(unlimited_storage&&) = default;
0430   unlimited_storage& operator=(unlimited_storage&&) = default;
0431 
0432   // TODO
0433   // template <class Allocator>
0434   // unlimited_storage(const unlimited_storage<Allocator>& s)
0435 
0436   template <class Iterable, class = detail::requires_iterable<Iterable>>
0437   explicit unlimited_storage(const Iterable& s) {
0438     using std::begin;
0439     using std::end;
0440     auto s_begin = begin(s);
0441     auto s_end = end(s);
0442     using V = typename std::iterator_traits<decltype(begin(s))>::value_type;
0443     // must be non-const to avoid msvc warning about if constexpr
0444     auto ti = buffer_type::template type_index<V>();
0445     auto nt = mp11::mp_size<typename buffer_type::types>::value;
0446     const std::size_t size = static_cast<std::size_t>(std::distance(s_begin, s_end));
0447     if (ti < nt)
0448       buffer_.template make<V>(size, s_begin);
0449     else
0450       buffer_.template make<double>(size, s_begin);
0451   }
0452 
0453   template <class Iterable, class = detail::requires_iterable<Iterable>>
0454   unlimited_storage& operator=(const Iterable& s) {
0455     *this = unlimited_storage(s);
0456     return *this;
0457   }
0458 
0459   allocator_type get_allocator() const { return buffer_.alloc; }
0460 
0461   void reset(std::size_t n) { buffer_.template make<U8>(n); }
0462 
0463   std::size_t size() const noexcept { return buffer_.size; }
0464 
0465   reference operator[](std::size_t i) noexcept { return {buffer_, i}; }
0466   const_reference operator[](std::size_t i) const noexcept { return {buffer_, i}; }
0467 
0468   bool operator==(const unlimited_storage& x) const noexcept {
0469     if (size() != x.size()) return false;
0470     return buffer_.visit([&x](const auto* p) {
0471       return x.buffer_.visit([p, n = x.size()](const auto* xp) {
0472         return std::equal(p, p + n, xp, detail::safe_equal{});
0473       });
0474     });
0475   }
0476 
0477   template <class Iterable>
0478   bool operator==(const Iterable& iterable) const {
0479     if (size() != iterable.size()) return false;
0480     return buffer_.visit([&iterable](const auto* p) {
0481       return std::equal(p, p + iterable.size(), std::begin(iterable),
0482                         detail::safe_equal{});
0483     });
0484   }
0485 
0486   unlimited_storage& operator*=(const double x) {
0487     buffer_.visit(multiplier(), buffer_, x);
0488     return *this;
0489   }
0490 
0491   iterator begin() noexcept { return {&buffer_, 0}; }
0492   iterator end() noexcept { return {&buffer_, size()}; }
0493   const_iterator begin() const noexcept { return {&buffer_, 0}; }
0494   const_iterator end() const noexcept { return {&buffer_, size()}; }
0495 
0496   /// implementation detail; used by unit tests, not part of generic storage interface
0497   template <class T>
0498   unlimited_storage(std::size_t s, const T* p, const allocator_type& a = {})
0499       : buffer_(std::move(a)) {
0500     buffer_.template make<T>(s, p);
0501   }
0502 
0503   template <class Archive>
0504   void serialize(Archive& ar, unsigned /* version */) {
0505     if (Archive::is_loading::value) {
0506       buffer_type tmp(buffer_.alloc);
0507       std::size_t size;
0508       ar& make_nvp("type", tmp.type);
0509       ar& make_nvp("size", size);
0510       tmp.visit([this, size](auto* tp) {
0511         assert(tp == nullptr);
0512         using T = std::decay_t<decltype(*tp)>;
0513         buffer_.template make<T>(size);
0514       });
0515     } else {
0516       ar& make_nvp("type", buffer_.type);
0517       ar& make_nvp("size", buffer_.size);
0518     }
0519     buffer_.visit([this, &ar](auto* tp) {
0520       auto w = detail::make_array_wrapper(tp, this->buffer_.size);
0521       ar& make_nvp("buffer", w);
0522     });
0523   }
0524 
0525 private:
0526   struct incrementor {
0527     template <class T>
0528     void operator()(T* tp, buffer_type& b, std::size_t i) {
0529       assert(tp && i < b.size);
0530       if (!detail::safe_increment(tp[i])) {
0531         using U = detail::next_type<typename buffer_type::types, T>;
0532         b.template make<U>(b.size, tp);
0533         ++static_cast<U*>(b.ptr)[i];
0534       }
0535     }
0536 
0537     void operator()(large_int* tp, buffer_type&, std::size_t i) { ++tp[i]; }
0538 
0539     void operator()(double* tp, buffer_type&, std::size_t i) { ++tp[i]; }
0540   };
0541 
0542   struct adder {
0543     template <class U>
0544     void operator()(double* tp, buffer_type&, std::size_t i, const U& x) {
0545       tp[i] += static_cast<double>(x);
0546     }
0547 
0548     void operator()(large_int* tp, buffer_type&, std::size_t i, const large_int& x) {
0549       tp[i] += x; // potentially adding large_int to itself is safe
0550     }
0551 
0552     template <class T, class U>
0553     void operator()(T* tp, buffer_type& b, std::size_t i, const U& x) {
0554       is_x_integral(std::is_integral<U>{}, tp, b, i, x);
0555     }
0556 
0557     template <class T, class U>
0558     void is_x_integral(std::false_type, T* tp, buffer_type& b, std::size_t i,
0559                        const U& x) {
0560       // x could be reference to buffer we manipulate, make copy before changing buffer
0561       const auto v = static_cast<double>(x);
0562       b.template make<double>(b.size, tp);
0563       operator()(static_cast<double*>(b.ptr), b, i, v);
0564     }
0565 
0566     template <class T>
0567     void is_x_integral(std::false_type, T* tp, buffer_type& b, std::size_t i,
0568                        const large_int& x) {
0569       // x could be reference to buffer we manipulate, make copy before changing buffer
0570       const auto v = static_cast<large_int>(x);
0571       b.template make<large_int>(b.size, tp);
0572       operator()(static_cast<large_int*>(b.ptr), b, i, v);
0573     }
0574 
0575     template <class T, class U>
0576     void is_x_integral(std::true_type, T* tp, buffer_type& b, std::size_t i, const U& x) {
0577       is_x_unsigned(std::is_unsigned<U>{}, tp, b, i, x);
0578     }
0579 
0580     template <class T, class U>
0581     void is_x_unsigned(std::false_type, T* tp, buffer_type& b, std::size_t i,
0582                        const U& x) {
0583       if (x >= 0)
0584         is_x_unsigned(std::true_type{}, tp, b, i, detail::make_unsigned(x));
0585       else
0586         is_x_integral(std::false_type{}, tp, b, i, static_cast<double>(x));
0587     }
0588 
0589     template <class T, class U>
0590     void is_x_unsigned(std::true_type, T* tp, buffer_type& b, std::size_t i, const U& x) {
0591       if (detail::safe_radd(tp[i], x)) return;
0592       // x could be reference to buffer we manipulate, need to convert to value
0593       const auto y = x;
0594       using TN = detail::next_type<typename buffer_type::types, T>;
0595       b.template make<TN>(b.size, tp);
0596       is_x_unsigned(std::true_type{}, static_cast<TN*>(b.ptr), b, i, y);
0597     }
0598 
0599     template <class U>
0600     void is_x_unsigned(std::true_type, large_int* tp, buffer_type&, std::size_t i,
0601                        const U& x) {
0602       tp[i] += x;
0603     }
0604   };
0605 
0606   struct multiplier {
0607     template <class T>
0608     void operator()(T* tp, buffer_type& b, const double x) {
0609       // potential lossy conversion that cannot be avoided
0610       b.template make<double>(b.size, tp);
0611       operator()(static_cast<double*>(b.ptr), b, x);
0612     }
0613 
0614     void operator()(double* tp, buffer_type& b, const double x) {
0615       for (auto end = tp + b.size; tp != end; ++tp) *tp *= x;
0616     }
0617 
0618     template <class T>
0619     void operator()(T* tp, buffer_type& b, std::size_t i, const double x) {
0620       b.template make<double>(b.size, tp);
0621       operator()(static_cast<double*>(b.ptr), b, i, x);
0622     }
0623 
0624     void operator()(double* tp, buffer_type&, std::size_t i, const double x) {
0625       tp[i] *= static_cast<double>(x);
0626     }
0627   };
0628 
0629   mutable buffer_type buffer_;
0630   friend struct unsafe_access;
0631 };
0632 
0633 } // namespace histogram
0634 } // namespace boost
0635 
0636 #endif