boost/histogram/detail/large_int.hpp

   1 // Copyright 2018-2019 Hans Dembinski
   2 //
   3 // Distributed under the Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_HISTOGRAM_DETAIL_LARGE_INT_HPP
   8 #define BOOST_HISTOGRAM_DETAIL_LARGE_INT_HPP
   9
  10 #include <boost/assert.hpp>
  11 #include <boost/histogram/detail/operators.hpp>
  12 #include <boost/histogram/detail/safe_comparison.hpp>
  13 #include <boost/mp11/algorithm.hpp>
  14 #include <boost/mp11/function.hpp>
  15 #include <boost/mp11/list.hpp>
  16 #include <boost/mp11/utility.hpp>
  17 #include <cmath>
  18 #include <cstdint>
  19 #include <limits>
  20 #include <type_traits>
  21 #include <utility>
  22 #include <vector>
  23
  24 namespace boost {
  25 namespace histogram {
  26 namespace detail {
  27
  28 template <class T>
  29 using is_unsigned_integral = mp11::mp_and<std::is_integral<T>, std::is_unsigned<T>>;
  30
  31 template <class T>
  32 bool safe_increment(T& t) {
  33   if (t < (std::numeric_limits<T>::max)()) {
  34     ++t;
  35     return true;
  36   }
  37   return false;
  38 }
  39
  40 template <class T, class U>
  41 bool safe_radd(T& t, const U& u) {
  42   static_assert(is_unsigned_integral<T>::value, "T must be unsigned integral type");
  43   static_assert(is_unsigned_integral<U>::value, "T must be unsigned integral type");
  44   if (static_cast<T>((std::numeric_limits<T>::max)() - t) >= u) {
  45     t += static_cast<T>(u); // static_cast to suppress conversion warning
  46     return true;
  47   }
  48   return false;
  49 }
  50
  51 // An integer type which can grow arbitrarily large (until memory is exhausted).
  52 // Use boost.multiprecision.cpp_int in your own code, it is much more sophisticated.
  53 // We use it only to reduce coupling between boost libs.
  54 template <class Allocator>
  55 struct large_int : totally_ordered<large_int<Allocator>, large_int<Allocator>>,
  56                    partially_ordered<large_int<Allocator>, void> {
  57   explicit large_int(const Allocator& a = {}) : data(1, 0, a) {}
  58   explicit large_int(std::uint64_t v, const Allocator& a = {}) : data(1, v, a) {}
  59
  60   large_int(const large_int&) = default;
  61   large_int& operator=(const large_int&) = default;
  62   large_int(large_int&&) = default;
  63   large_int& operator=(large_int&&) = default;
  64
  65   large_int& operator=(std::uint64_t o) {
  66     data = decltype(data)(1, o);
  67     return *this;
  68   }
  69
  70   large_int& operator++() {
  71     BOOST_ASSERT(data.size() > 0u);
  72     std::size_t i = 0;
  73     while (!safe_increment(data[i])) {
  74       data[i] = 0;
  75       ++i;
  76       if (i == data.size()) {
  77         data.push_back(1);
  78         break;
  79       }
  80     }
  81     return *this;
  82   }
  83
  84   large_int& operator+=(const large_int& o) {
  85     if (this == &o) {
  86       auto tmp{o};
  87       return operator+=(tmp);
  88     }
  89     bool carry = false;
  90     std::size_t i = 0;
  91     for (std::uint64_t oi : o.data) {
  92       auto& di = maybe_extend(i);
  93       if (carry) {
  94         if (safe_increment(oi))
  95           carry = false;
  96         else {
  97           ++i;
  98           continue;
  99         }
 100       }
 101       if (!safe_radd(di, oi)) {
 102         add_remainder(di, oi);
 103         carry = true;
 104       }
 105       ++i;
 106     }
 107     while (carry) {
 108       auto& di = maybe_extend(i);
 109       if (safe_increment(di)) break;
 110       di = 0;
 111       ++i;
 112     }
 113     return *this;
 114   }
 115
 116   large_int& operator+=(std::uint64_t o) {
 117     BOOST_ASSERT(data.size() > 0u);
 118     if (safe_radd(data[0], o)) return *this;
 119     add_remainder(data[0], o);
 120     // carry the one, data may grow several times
 121     std::size_t i = 1;
 122     while (true) {
 123       auto& di = maybe_extend(i);
 124       if (safe_increment(di)) break;
 125       di = 0;
 126       ++i;
 127     }
 128     return *this;
 129   }
 130
 131   explicit operator double() const noexcept {
 132     BOOST_ASSERT(data.size() > 0u);
 133     double result = static_cast<double>(data[0]);
 134     std::size_t i = 0;
 135     while (++i < data.size())
 136       result += static_cast<double>(data[i]) * std::pow(2.0, i * 64);
 137     return result;
 138   }
 139
 140   bool operator<(const large_int& o) const noexcept {
 141     BOOST_ASSERT(data.size() > 0u);
 142     BOOST_ASSERT(o.data.size() > 0u);
 143     // no leading zeros allowed
 144     BOOST_ASSERT(data.size() == 1 || data.back() > 0u);
 145     BOOST_ASSERT(o.data.size() == 1 || o.data.back() > 0u);
 146     if (data.size() < o.data.size()) return true;
 147     if (data.size() > o.data.size()) return false;
 148     auto s = data.size();
 149     while (s > 0u) {
 150       --s;
 151       if (data[s] < o.data[s]) return true;
 152       if (data[s] > o.data[s]) return false;
 153     }
 154     return false; // args are equal
 155   }
 156
 157   bool operator==(const large_int& o) const noexcept {
 158     BOOST_ASSERT(data.size() > 0u);
 159     BOOST_ASSERT(o.data.size() > 0u);
 160     // no leading zeros allowed
 161     BOOST_ASSERT(data.size() == 1 || data.back() > 0u);
 162     BOOST_ASSERT(o.data.size() == 1 || o.data.back() > 0u);
 163     if (data.size() != o.data.size()) return false;
 164     return std::equal(data.begin(), data.end(), o.data.begin());
 165   }
 166
 167   template <class U>
 168   std::enable_if_t<std::is_integral<U>::value, bool> operator<(const U& o) const
 169       noexcept {
 170     BOOST_ASSERT(data.size() > 0u);
 171     return data.size() == 1 && safe_less()(data[0], o);
 172   }
 173
 174   template <class U>
 175   std::enable_if_t<std::is_integral<U>::value, bool> operator>(const U& o) const
 176       noexcept {
 177     BOOST_ASSERT(data.size() > 0u);
 178     BOOST_ASSERT(data.size() == 1 || data.back() > 0u); // no leading zeros allowed
 179     return data.size() > 1 || safe_less()(o, data[0]);
 180   }
 181
 182   template <class U>
 183   std::enable_if_t<std::is_integral<U>::value, bool> operator==(const U& o) const
 184       noexcept {
 185     BOOST_ASSERT(data.size() > 0u);
 186     return data.size() == 1 && safe_equal()(data[0], o);
 187   }
 188
 189   template <class U>
 190   std::enable_if_t<std::is_floating_point<U>::value, bool> operator<(const U& o) const
 191       noexcept {
 192     return operator double() < o;
 193   }
 194   template <class U>
 195   std::enable_if_t<std::is_floating_point<U>::value, bool> operator>(const U& o) const
 196       noexcept {
 197     return operator double() > o;
 198   }
 199   template <class U>
 200   std::enable_if_t<std::is_floating_point<U>::value, bool> operator==(const U& o) const
 201       noexcept {
 202     return operator double() == o;
 203   }
 204
 205   std::uint64_t& maybe_extend(std::size_t i) {
 206     while (i >= data.size()) data.push_back(0);
 207     return data[i];
 208   }
 209
 210   static void add_remainder(std::uint64_t& d, const std::uint64_t o) noexcept {
 211     BOOST_ASSERT(d > 0u);
 212     // in decimal system it would look like this:
 213     // 8 + 8 = 6 = 8 - (9 - 8) - 1
 214     // 9 + 1 = 0 = 9 - (9 - 1) - 1
 215     auto tmp = (std::numeric_limits<std::uint64_t>::max)();
 216     tmp -= o;
 217     --d -= tmp;
 218   }
 219
 220   template <class Archive>
 221   void serialize(Archive& ar, unsigned /* version */) {
 222     ar& make_nvp("data", data);
 223   }
 224
 225   std::vector<std::uint64_t, Allocator> data;
 226 };
 227
 228 } // namespace detail
 229 } // namespace histogram
 230 } // namespace boost
 231
 232 #endif