histogram/detail/large_int.hpp

0001 // Copyright 2018-2019 Hans Dembinski
0002 //
0003 // Distributed under the Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_HISTOGRAM_DETAIL_LARGE_INT_HPP
0008 #define BOOST_HISTOGRAM_DETAIL_LARGE_INT_HPP
0009
0010 #include <boost/histogram/detail/operators.hpp>
0011 #include <boost/histogram/detail/safe_comparison.hpp>
0012 #include <boost/mp11/algorithm.hpp>
0013 #include <boost/mp11/function.hpp>
0014 #include <boost/mp11/list.hpp>
0015 #include <boost/mp11/utility.hpp>
0016 #include <cassert>
0017 #include <cmath>
0018 #include <cstdint>
0019 #include <limits>
0020 #include <type_traits>
0021 #include <utility>
0022 #include <vector>
0023
0024 namespace boost {
0025 namespace histogram {
0026 namespace detail {
0027
0028 template <class T>
0029 using is_unsigned_integral = mp11::mp_and<std::is_integral<T>, std::is_unsigned<T>>;
0030
0031 template <class T>
0032 bool safe_increment(T& t) {
0033   if (t < (std::numeric_limits<T>::max)()) {
0034     ++t;
0035     return true;
0036   }
0037   return false;
0038 }
0039
0040 template <class T, class U>
0041 bool safe_radd(T& t, const U& u) {
0042   static_assert(is_unsigned_integral<T>::value, "T must be unsigned integral type");
0043   static_assert(is_unsigned_integral<U>::value, "T must be unsigned integral type");
0044   if (static_cast<T>((std::numeric_limits<T>::max)() - t) >= u) {
0045     t += static_cast<T>(u); // static_cast to suppress conversion warning
0046     return true;
0047   }
0048   return false;
0049 }
0050
0051 // An integer type which can grow arbitrarily large (until memory is exhausted).
0052 // Use boost.multiprecision.cpp_int in your own code, it is much more sophisticated.
0053 // We use it only to reduce coupling between boost libs.
0054 template <class Allocator>
0055 struct large_int : totally_ordered<large_int<Allocator>, large_int<Allocator>>,
0056                    partially_ordered<large_int<Allocator>, void> {
0057   explicit large_int(const Allocator& a = {}) : data(1, 0, a) {}
0058   explicit large_int(std::uint64_t v, const Allocator& a = {}) : data(1, v, a) {}
0059
0060   large_int(const large_int&) = default;
0061   large_int& operator=(const large_int&) = default;
0062   large_int(large_int&&) = default;
0063   large_int& operator=(large_int&&) = default;
0064
0065   large_int& operator=(std::uint64_t o) {
0066     data = decltype(data)(1, o);
0067     return *this;
0068   }
0069
0070   large_int& operator++() {
0071     assert(data.size() > 0u);
0072     std::size_t i = 0;
0073     while (!safe_increment(data[i])) {
0074       data[i] = 0;
0075       ++i;
0076       if (i == data.size()) {
0077         data.push_back(1);
0078         break;
0079       }
0080     }
0081     return *this;
0082   }
0083
0084   large_int& operator+=(const large_int& o) {
0085     if (this == &o) {
0086       auto tmp{o};
0087       return operator+=(tmp);
0088     }
0089     bool carry = false;
0090     std::size_t i = 0;
0091     for (std::uint64_t oi : o.data) {
0092       auto& di = maybe_extend(i);
0093       if (carry) {
0094         if (safe_increment(oi))
0095           carry = false;
0096         else {
0097           ++i;
0098           continue;
0099         }
0100       }
0101       if (!safe_radd(di, oi)) {
0102         add_remainder(di, oi);
0103         carry = true;
0104       }
0105       ++i;
0106     }
0107     while (carry) {
0108       auto& di = maybe_extend(i);
0109       if (safe_increment(di)) break;
0110       di = 0;
0111       ++i;
0112     }
0113     return *this;
0114   }
0115
0116   large_int& operator+=(std::uint64_t o) {
0117     assert(data.size() > 0u);
0118     if (safe_radd(data[0], o)) return *this;
0119     add_remainder(data[0], o);
0120     // carry the one, data may grow several times
0121     std::size_t i = 1;
0122     while (true) {
0123       auto& di = maybe_extend(i);
0124       if (safe_increment(di)) break;
0125       di = 0;
0126       ++i;
0127     }
0128     return *this;
0129   }
0130
0131   explicit operator double() const noexcept {
0132     assert(data.size() > 0u);
0133     double result = static_cast<double>(data[0]);
0134     std::size_t i = 0;
0135     while (++i < data.size())
0136       result += static_cast<double>(data[i]) * std::pow(2.0, i * 64);
0137     return result;
0138   }
0139
0140   bool operator<(const large_int& o) const noexcept {
0141     assert(data.size() > 0u);
0142     assert(o.data.size() > 0u);
0143     // no leading zeros allowed
0144     assert(data.size() == 1 || data.back() > 0u);
0145     assert(o.data.size() == 1 || o.data.back() > 0u);
0146     if (data.size() < o.data.size()) return true;
0147     if (data.size() > o.data.size()) return false;
0148     auto s = data.size();
0149     while (s > 0u) {
0150       --s;
0151       if (data[s] < o.data[s]) return true;
0152       if (data[s] > o.data[s]) return false;
0153     }
0154     return false; // args are equal
0155   }
0156
0157   bool operator==(const large_int& o) const noexcept {
0158     assert(data.size() > 0u);
0159     assert(o.data.size() > 0u);
0160     // no leading zeros allowed
0161     assert(data.size() == 1 || data.back() > 0u);
0162     assert(o.data.size() == 1 || o.data.back() > 0u);
0163     if (data.size() != o.data.size()) return false;
0164     return std::equal(data.begin(), data.end(), o.data.begin());
0165   }
0166
0167   template <class U>
0168   std::enable_if_t<std::is_integral<U>::value, bool> operator<(
0169       const U& o) const noexcept {
0170     assert(data.size() > 0u);
0171     return data.size() == 1 && safe_less()(data[0], o);
0172   }
0173
0174   template <class U>
0175   std::enable_if_t<std::is_integral<U>::value, bool> operator>(
0176       const U& o) const noexcept {
0177     assert(data.size() > 0u);
0178     assert(data.size() == 1 || data.back() > 0u); // no leading zeros allowed
0179     return data.size() > 1 || safe_less()(o, data[0]);
0180   }
0181
0182   template <class U>
0183   std::enable_if_t<std::is_integral<U>::value, bool> operator==(
0184       const U& o) const noexcept {
0185     assert(data.size() > 0u);
0186     return data.size() == 1 && safe_equal()(data[0], o);
0187   }
0188
0189   template <class U>
0190   std::enable_if_t<std::is_floating_point<U>::value, bool> operator<(
0191       const U& o) const noexcept {
0192     return operator double() < o;
0193   }
0194
0195   template <class U>
0196   std::enable_if_t<std::is_floating_point<U>::value, bool> operator>(
0197       const U& o) const noexcept {
0198     return operator double() > o;
0199   }
0200
0201   template <class U>
0202   std::enable_if_t<std::is_floating_point<U>::value, bool> operator==(
0203       const U& o) const noexcept {
0204     return operator double() == o;
0205   }
0206
0207   template <class U>
0208   std::enable_if_t<
0209       (!std::is_arithmetic<U>::value && std::is_convertible<U, double>::value), bool>
0210   operator<(const U& o) const noexcept {
0211     return operator double() < o;
0212   }
0213
0214   template <class U>
0215   std::enable_if_t<
0216       (!std::is_arithmetic<U>::value && std::is_convertible<U, double>::value), bool>
0217   operator>(const U& o) const noexcept {
0218     return operator double() > o;
0219   }
0220
0221   template <class U>
0222   std::enable_if_t<
0223       (!std::is_arithmetic<U>::value && std::is_convertible<U, double>::value), bool>
0224   operator==(const U& o) const noexcept {
0225     return operator double() == o;
0226   }
0227
0228   std::uint64_t& maybe_extend(std::size_t i) {
0229     while (i >= data.size()) data.push_back(0);
0230     return data[i];
0231   }
0232
0233   static void add_remainder(std::uint64_t& d, const std::uint64_t o) noexcept {
0234     assert(d > 0u);
0235     // in decimal system it would look like this:
0236     // 8 + 8 = 6 = 8 - (9 - 8) - 1
0237     // 9 + 1 = 0 = 9 - (9 - 1) - 1
0238     auto tmp = (std::numeric_limits<std::uint64_t>::max)();
0239     tmp -= o;
0240     --d -= tmp;
0241   }
0242
0243   template <class Archive>
0244   void serialize(Archive& ar, unsigned /* version */) {
0245     ar& make_nvp("data", data);
0246   }
0247
0248   std::vector<std::uint64_t, Allocator> data;
0249 };
0250
0251 } // namespace detail
0252 } // namespace histogram
0253 } // namespace boost
0254
0255 #endif