integer.hpp 8.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246
  1. ///////////////////////////////////////////////////////////////
  2. // Copyright 2012 John Maddock. Distributed under the Boost
  3. // Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
  5. #ifndef BOOST_MP_INTEGER_HPP
  6. #define BOOST_MP_INTEGER_HPP
  7. #include <boost/multiprecision/cpp_int.hpp>
  8. #include <boost/multiprecision/detail/bitscan.hpp>
  9. namespace boost {
  10. namespace multiprecision {
  11. template <class Integer, class I2>
  12. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && boost::multiprecision::detail::is_integral<I2>::value, Integer&>::type
  13. multiply(Integer& result, const I2& a, const I2& b)
  14. {
  15. return result = static_cast<Integer>(a) * static_cast<Integer>(b);
  16. }
  17. template <class Integer, class I2>
  18. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && boost::multiprecision::detail::is_integral<I2>::value, Integer&>::type
  19. add(Integer& result, const I2& a, const I2& b)
  20. {
  21. return result = static_cast<Integer>(a) + static_cast<Integer>(b);
  22. }
  23. template <class Integer, class I2>
  24. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && boost::multiprecision::detail::is_integral<I2>::value, Integer&>::type
  25. subtract(Integer& result, const I2& a, const I2& b)
  26. {
  27. return result = static_cast<Integer>(a) - static_cast<Integer>(b);
  28. }
  29. template <class Integer>
  30. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value>::type divide_qr(const Integer& x, const Integer& y, Integer& q, Integer& r)
  31. {
  32. q = x / y;
  33. r = x % y;
  34. }
  35. template <class I1, class I2>
  36. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<I1>::value && boost::multiprecision::detail::is_integral<I2>::value, I2>::type integer_modulus(const I1& x, I2 val)
  37. {
  38. return static_cast<I2>(x % val);
  39. }
  40. namespace detail {
  41. //
  42. // Figure out the kind of integer that has twice as many bits as some builtin
  43. // integer type I. Use a native type if we can (including types which may not
  44. // be recognised by boost::int_t because they're larger than boost::long_long_type),
  45. // otherwise synthesize a cpp_int to do the job.
  46. //
  47. template <class I>
  48. struct double_integer
  49. {
  50. static constexpr const unsigned int_t_digits =
  51. 2 * sizeof(I) <= sizeof(boost::long_long_type) ? std::numeric_limits<I>::digits * 2 : 1;
  52. using type = typename std::conditional<
  53. 2 * sizeof(I) <= sizeof(boost::long_long_type),
  54. typename std::conditional<
  55. boost::multiprecision::detail::is_signed<I>::value && boost::multiprecision::detail::is_integral<I>::value,
  56. typename boost::int_t<int_t_digits>::least,
  57. typename boost::uint_t<int_t_digits>::least>::type,
  58. typename std::conditional<
  59. 2 * sizeof(I) <= sizeof(double_limb_type),
  60. typename std::conditional<
  61. boost::multiprecision::detail::is_signed<I>::value && boost::multiprecision::detail::is_integral<I>::value,
  62. signed_double_limb_type,
  63. double_limb_type>::type,
  64. number<cpp_int_backend<sizeof(I) * CHAR_BIT * 2, sizeof(I) * CHAR_BIT * 2, (boost::multiprecision::detail::is_signed<I>::value ? signed_magnitude : unsigned_magnitude), unchecked, void> > >::type>::type;
  65. };
  66. } // namespace detail
  67. template <class I1, class I2, class I3>
  68. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<I1>::value && boost::multiprecision::detail::is_unsigned<I2>::value && is_integral<I3>::value, I1>::type
  69. powm(const I1& a, I2 b, I3 c)
  70. {
  71. using double_type = typename detail::double_integer<I1>::type;
  72. I1 x(1), y(a);
  73. double_type result(0);
  74. while (b > 0)
  75. {
  76. if (b & 1)
  77. {
  78. multiply(result, x, y);
  79. x = integer_modulus(result, c);
  80. }
  81. multiply(result, y, y);
  82. y = integer_modulus(result, c);
  83. b >>= 1;
  84. }
  85. return x % c;
  86. }
  87. template <class I1, class I2, class I3>
  88. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<I1>::value && boost::multiprecision::detail::is_signed<I2>::value && boost::multiprecision::detail::is_integral<I2>::value && boost::multiprecision::detail::is_integral<I3>::value, I1>::type
  89. powm(const I1& a, I2 b, I3 c)
  90. {
  91. if (b < 0)
  92. {
  93. BOOST_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
  94. }
  95. return powm(a, static_cast<typename boost::multiprecision::detail::make_unsigned<I2>::type>(b), c);
  96. }
  97. template <class Integer>
  98. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, unsigned>::type lsb(const Integer& val)
  99. {
  100. if (val <= 0)
  101. {
  102. if (val == 0)
  103. {
  104. BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
  105. }
  106. else
  107. {
  108. BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
  109. }
  110. }
  111. return detail::find_lsb(val);
  112. }
  113. template <class Integer>
  114. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, unsigned>::type msb(Integer val)
  115. {
  116. if (val <= 0)
  117. {
  118. if (val == 0)
  119. {
  120. BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
  121. }
  122. else
  123. {
  124. BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
  125. }
  126. }
  127. return detail::find_msb(val);
  128. }
  129. template <class Integer>
  130. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, bool>::type bit_test(const Integer& val, unsigned index)
  131. {
  132. Integer mask = 1;
  133. if (index >= sizeof(Integer) * CHAR_BIT)
  134. return 0;
  135. if (index)
  136. mask <<= index;
  137. return val & mask ? true : false;
  138. }
  139. template <class Integer>
  140. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer&>::type bit_set(Integer& val, unsigned index)
  141. {
  142. Integer mask = 1;
  143. if (index >= sizeof(Integer) * CHAR_BIT)
  144. return val;
  145. if (index)
  146. mask <<= index;
  147. val |= mask;
  148. return val;
  149. }
  150. template <class Integer>
  151. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer&>::type bit_unset(Integer& val, unsigned index)
  152. {
  153. Integer mask = 1;
  154. if (index >= sizeof(Integer) * CHAR_BIT)
  155. return val;
  156. if (index)
  157. mask <<= index;
  158. val &= ~mask;
  159. return val;
  160. }
  161. template <class Integer>
  162. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer&>::type bit_flip(Integer& val, unsigned index)
  163. {
  164. Integer mask = 1;
  165. if (index >= sizeof(Integer) * CHAR_BIT)
  166. return val;
  167. if (index)
  168. mask <<= index;
  169. val ^= mask;
  170. return val;
  171. }
  172. template <class Integer>
  173. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer>::type sqrt(const Integer& x, Integer& r)
  174. {
  175. //
  176. // This is slow bit-by-bit integer square root, see for example
  177. // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
  178. // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
  179. // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
  180. // at some point.
  181. //
  182. Integer s = 0;
  183. if (x == 0)
  184. {
  185. r = 0;
  186. return s;
  187. }
  188. int g = msb(x);
  189. if (g == 0)
  190. {
  191. r = 1;
  192. return s;
  193. }
  194. Integer t = 0;
  195. r = x;
  196. g /= 2;
  197. bit_set(s, g);
  198. bit_set(t, 2 * g);
  199. r = x - t;
  200. --g;
  201. do
  202. {
  203. t = s;
  204. t <<= g + 1;
  205. bit_set(t, 2 * g);
  206. if (t <= r)
  207. {
  208. bit_set(s, g);
  209. r -= t;
  210. }
  211. --g;
  212. } while (g >= 0);
  213. return s;
  214. }
  215. template <class Integer>
  216. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value, Integer>::type sqrt(const Integer& x)
  217. {
  218. Integer r(0);
  219. return sqrt(x, r);
  220. }
  221. }} // namespace boost::multiprecision
  222. #endif