integer.hpp 6.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240
  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 http://www.boost.org/LICENSE_1_
  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. typename enable_if_c<is_integral<Integer>::value && 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. typename enable_if_c<is_integral<Integer>::value && 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. typename enable_if_c<is_integral<Integer>::value && 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. typename enable_if_c<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. typename enable_if_c<is_integral<I1>::value && 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 long long),
  45. // otherwise synthesize a cpp_int to do the job.
  46. //
  47. template <class I>
  48. struct double_integer
  49. {
  50. static const unsigned int_t_digits =
  51. 2 * sizeof(I) <= sizeof(long long) ? std::numeric_limits<I>::digits * 2 : 1;
  52. typedef typename mpl::if_c<
  53. 2 * sizeof(I) <= sizeof(long long),
  54. typename mpl::if_c<
  55. is_signed<I>::value,
  56. typename boost::int_t<int_t_digits>::least,
  57. typename boost::uint_t<int_t_digits>::least
  58. >::type,
  59. typename mpl::if_c<
  60. 2 * sizeof(I) <= sizeof(double_limb_type),
  61. typename mpl::if_c<
  62. is_signed<I>::value,
  63. signed_double_limb_type,
  64. double_limb_type
  65. >::type,
  66. number<cpp_int_backend<sizeof(I)*CHAR_BIT*2, sizeof(I)*CHAR_BIT*2, (is_signed<I>::value ? signed_magnitude : unsigned_magnitude), unchecked, void> >
  67. >::type
  68. >::type type;
  69. };
  70. }
  71. template <class I1, class I2, class I3>
  72. typename enable_if_c<is_integral<I1>::value && is_integral<I2>::value && is_integral<I3>::value, I1>::type
  73. powm(const I1& a, I2 b, I3 c)
  74. {
  75. typedef typename detail::double_integer<I1>::type double_type;
  76. I1 x(1), y(a);
  77. double_type result;
  78. while(b > 0)
  79. {
  80. if(b & 1)
  81. {
  82. multiply(result, x, y);
  83. x = integer_modulus(result, c);
  84. }
  85. multiply(result, y, y);
  86. y = integer_modulus(result, c);
  87. b >>= 1;
  88. }
  89. return x % c;
  90. }
  91. template <class Integer>
  92. typename enable_if_c<is_integral<Integer>::value, unsigned>::type lsb(const Integer& val)
  93. {
  94. if(val <= 0)
  95. {
  96. if(val == 0)
  97. {
  98. BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
  99. }
  100. else
  101. {
  102. BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
  103. }
  104. }
  105. return detail::find_lsb(val);
  106. }
  107. template <class Integer>
  108. typename enable_if_c<is_integral<Integer>::value, unsigned>::type msb(Integer val)
  109. {
  110. if(val <= 0)
  111. {
  112. if(val == 0)
  113. {
  114. BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
  115. }
  116. else
  117. {
  118. BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
  119. }
  120. }
  121. return detail::find_msb(val);
  122. }
  123. template <class Integer>
  124. typename enable_if_c<is_integral<Integer>::value, bool>::type bit_test(const Integer& val, unsigned index)
  125. {
  126. Integer mask = 1;
  127. if(index >= sizeof(Integer) * CHAR_BIT)
  128. return 0;
  129. if(index)
  130. mask <<= index;
  131. return val & mask ? true : false;
  132. }
  133. template <class Integer>
  134. typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_set(Integer& val, unsigned index)
  135. {
  136. Integer mask = 1;
  137. if(index >= sizeof(Integer) * CHAR_BIT)
  138. return val;
  139. if(index)
  140. mask <<= index;
  141. val |= mask;
  142. return val;
  143. }
  144. template <class Integer>
  145. typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_unset(Integer& val, unsigned index)
  146. {
  147. Integer mask = 1;
  148. if(index >= sizeof(Integer) * CHAR_BIT)
  149. return val;
  150. if(index)
  151. mask <<= index;
  152. val &= ~mask;
  153. return val;
  154. }
  155. template <class Integer>
  156. typename enable_if_c<is_integral<Integer>::value, Integer&>::type bit_flip(Integer& val, unsigned index)
  157. {
  158. Integer mask = 1;
  159. if(index >= sizeof(Integer) * CHAR_BIT)
  160. return val;
  161. if(index)
  162. mask <<= index;
  163. val ^= mask;
  164. return val;
  165. }
  166. template <class Integer>
  167. typename enable_if_c<is_integral<Integer>::value, Integer>::type sqrt(const Integer& x, Integer& r)
  168. {
  169. //
  170. // This is slow bit-by-bit integer square root, see for example
  171. // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
  172. // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
  173. // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
  174. // at some point.
  175. //
  176. Integer s = 0;
  177. if(x == 0)
  178. {
  179. r = 0;
  180. return s;
  181. }
  182. int g = msb(x);
  183. if(g == 0)
  184. {
  185. r = 1;
  186. return s;
  187. }
  188. Integer t;
  189. r = x;
  190. g /= 2;
  191. bit_set(s, g);
  192. bit_set(t, 2 * g);
  193. r = x - t;
  194. --g;
  195. do
  196. {
  197. t = s;
  198. t <<= g + 1;
  199. bit_set(t, 2 * g);
  200. if(t <= r)
  201. {
  202. bit_set(s, g);
  203. r -= t;
  204. }
  205. --g;
  206. }
  207. while(g >= 0);
  208. return s;
  209. }
  210. template <class Integer>
  211. typename enable_if_c<is_integral<Integer>::value, Integer>::type sqrt(const Integer& x)
  212. {
  213. Integer r;
  214. return sqrt(x, r);
  215. }
  216. }} // namespaces
  217. #endif