integer_ops.hpp 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474
  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_DETAIL_INTEGER_OPS_HPP
  6. #define BOOST_MP_DETAIL_INTEGER_OPS_HPP
  7. #include <boost/multiprecision/number.hpp>
  8. #include <boost/multiprecision/detail/no_exceptions_support.hpp>
  9. namespace boost { namespace multiprecision {
  10. namespace default_ops {
  11. template <class Backend>
  12. inline BOOST_MP_CXX14_CONSTEXPR void eval_qr(const Backend& x, const Backend& y, Backend& q, Backend& r)
  13. {
  14. eval_divide(q, x, y);
  15. eval_modulus(r, x, y);
  16. }
  17. template <class Backend, class Integer>
  18. inline BOOST_MP_CXX14_CONSTEXPR Integer eval_integer_modulus(const Backend& x, Integer val)
  19. {
  20. BOOST_MP_USING_ABS
  21. using default_ops::eval_convert_to;
  22. using default_ops::eval_modulus;
  23. using int_type = typename boost::multiprecision::detail::canonical<Integer, Backend>::type;
  24. Backend t;
  25. eval_modulus(t, x, static_cast<int_type>(val));
  26. Integer result(0);
  27. eval_convert_to(&result, t);
  28. return abs(result);
  29. }
  30. template <class B>
  31. inline BOOST_MP_CXX14_CONSTEXPR void eval_gcd(B& result, const B& a, const B& b)
  32. {
  33. using default_ops::eval_get_sign;
  34. using default_ops::eval_is_zero;
  35. using default_ops::eval_lsb;
  36. std::ptrdiff_t shift(0);
  37. B u(a), v(b);
  38. int s = eval_get_sign(u);
  39. /* GCD(0,x) := x */
  40. if (s < 0)
  41. {
  42. u.negate();
  43. }
  44. else if (s == 0)
  45. {
  46. result = v;
  47. return;
  48. }
  49. s = eval_get_sign(v);
  50. if (s < 0)
  51. {
  52. v.negate();
  53. }
  54. else if (s == 0)
  55. {
  56. result = u;
  57. return;
  58. }
  59. /* Let shift := lg K, where K is the greatest power of 2
  60. dividing both u and v. */
  61. std::size_t us = eval_lsb(u);
  62. std::size_t vs = eval_lsb(v);
  63. shift = static_cast<std::ptrdiff_t>((std::min)(us, vs));
  64. eval_right_shift(u, us);
  65. eval_right_shift(v, vs);
  66. do
  67. {
  68. /* Now u and v are both odd, so diff(u, v) is even.
  69. Let u = min(u, v), v = diff(u, v)/2. */
  70. s = u.compare(v);
  71. if (s > 0)
  72. u.swap(v);
  73. if (s == 0)
  74. break;
  75. eval_subtract(v, u);
  76. vs = eval_lsb(v);
  77. eval_right_shift(v, vs);
  78. } while (true);
  79. result = u;
  80. eval_left_shift(result, shift);
  81. }
  82. template <class B>
  83. inline BOOST_MP_CXX14_CONSTEXPR void eval_lcm(B& result, const B& a, const B& b)
  84. {
  85. using ui_type = typename std::tuple_element<0, typename B::unsigned_types>::type;
  86. B t;
  87. eval_gcd(t, a, b);
  88. if (eval_is_zero(t))
  89. {
  90. result = static_cast<ui_type>(0);
  91. }
  92. else
  93. {
  94. eval_divide(result, a, t);
  95. eval_multiply(result, b);
  96. }
  97. if (eval_get_sign(result) < 0)
  98. result.negate();
  99. }
  100. } // namespace default_ops
  101. template <class Backend, expression_template_option ExpressionTemplates>
  102. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer>::type
  103. divide_qr(const number<Backend, ExpressionTemplates>& x, const number<Backend, ExpressionTemplates>& y,
  104. number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
  105. {
  106. using default_ops::eval_qr;
  107. eval_qr(x.backend(), y.backend(), q.backend(), r.backend());
  108. }
  109. template <class Backend, expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4>
  110. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer>::type
  111. divide_qr(const number<Backend, ExpressionTemplates>& x, const multiprecision::detail::expression<tag, A1, A2, A3, A4>& y,
  112. number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
  113. {
  114. divide_qr(x, number<Backend, ExpressionTemplates>(y), q, r);
  115. }
  116. template <class tag, class A1, class A2, class A3, class A4, class Backend, expression_template_option ExpressionTemplates>
  117. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer>::type
  118. divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const number<Backend, ExpressionTemplates>& y,
  119. number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
  120. {
  121. divide_qr(number<Backend, ExpressionTemplates>(x), y, q, r);
  122. }
  123. template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b, class Backend, expression_template_option ExpressionTemplates>
  124. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer>::type
  125. divide_qr(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, const multiprecision::detail::expression<tagb, A1b, A2b, A3b, A4b>& y,
  126. number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r)
  127. {
  128. divide_qr(number<Backend, ExpressionTemplates>(x), number<Backend, ExpressionTemplates>(y), q, r);
  129. }
  130. template <class Backend, expression_template_option ExpressionTemplates, class Integer>
  131. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && (number_category<Backend>::value == number_kind_integer), Integer>::type
  132. integer_modulus(const number<Backend, ExpressionTemplates>& x, Integer val)
  133. {
  134. using default_ops::eval_integer_modulus;
  135. return eval_integer_modulus(x.backend(), val);
  136. }
  137. template <class tag, class A1, class A2, class A3, class A4, class Integer>
  138. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<Integer>::value && (number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer), Integer>::type
  139. integer_modulus(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, Integer val)
  140. {
  141. using result_type = typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type;
  142. return integer_modulus(result_type(x), val);
  143. }
  144. template <class Backend, expression_template_option ExpressionTemplates>
  145. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, std::size_t>::type
  146. lsb(const number<Backend, ExpressionTemplates>& x)
  147. {
  148. using default_ops::eval_lsb;
  149. return eval_lsb(x.backend());
  150. }
  151. template <class tag, class A1, class A2, class A3, class A4>
  152. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, std::size_t>::type
  153. lsb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
  154. {
  155. using number_type = typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type;
  156. number_type n(x);
  157. using default_ops::eval_lsb;
  158. return eval_lsb(n.backend());
  159. }
  160. template <class Backend, expression_template_option ExpressionTemplates>
  161. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, std::size_t>::type
  162. msb(const number<Backend, ExpressionTemplates>& x)
  163. {
  164. using default_ops::eval_msb;
  165. return eval_msb(x.backend());
  166. }
  167. template <class tag, class A1, class A2, class A3, class A4>
  168. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, std::size_t>::type
  169. msb(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x)
  170. {
  171. using number_type = typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type;
  172. number_type n(x);
  173. using default_ops::eval_msb;
  174. return eval_msb(n.backend());
  175. }
  176. template <class Backend, expression_template_option ExpressionTemplates>
  177. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, bool>::type
  178. bit_test(const number<Backend, ExpressionTemplates>& x, std::size_t index)
  179. {
  180. using default_ops::eval_bit_test;
  181. return eval_bit_test(x.backend(), index);
  182. }
  183. template <class tag, class A1, class A2, class A3, class A4>
  184. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_integer, bool>::type
  185. bit_test(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& x, std::size_t index)
  186. {
  187. using number_type = typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type;
  188. number_type n(x);
  189. using default_ops::eval_bit_test;
  190. return eval_bit_test(n.backend(), index);
  191. }
  192. template <class Backend, expression_template_option ExpressionTemplates>
  193. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
  194. bit_set(number<Backend, ExpressionTemplates>& x, std::size_t index)
  195. {
  196. using default_ops::eval_bit_set;
  197. eval_bit_set(x.backend(), index);
  198. return x;
  199. }
  200. template <class Backend, expression_template_option ExpressionTemplates>
  201. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
  202. bit_unset(number<Backend, ExpressionTemplates>& x, std::size_t index)
  203. {
  204. using default_ops::eval_bit_unset;
  205. eval_bit_unset(x.backend(), index);
  206. return x;
  207. }
  208. template <class Backend, expression_template_option ExpressionTemplates>
  209. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<number_category<Backend>::value == number_kind_integer, number<Backend, ExpressionTemplates>&>::type
  210. bit_flip(number<Backend, ExpressionTemplates>& x, std::size_t index)
  211. {
  212. using default_ops::eval_bit_flip;
  213. eval_bit_flip(x.backend(), index);
  214. return x;
  215. }
  216. namespace default_ops {
  217. //
  218. // Within powm, we need a type with twice as many digits as the argument type, define
  219. // a traits class to obtain that type:
  220. //
  221. template <class Backend>
  222. struct double_precision_type
  223. {
  224. using type = Backend;
  225. };
  226. //
  227. // If the exponent is a signed integer type, then we need to
  228. // check the value is positive:
  229. //
  230. template <class Backend>
  231. inline BOOST_MP_CXX14_CONSTEXPR void check_sign_of_backend(const Backend& v, const std::integral_constant<bool, true>)
  232. {
  233. if (eval_get_sign(v) < 0)
  234. {
  235. BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
  236. }
  237. }
  238. template <class Backend>
  239. inline BOOST_MP_CXX14_CONSTEXPR void check_sign_of_backend(const Backend&, const std::integral_constant<bool, false>) {}
  240. //
  241. // Calculate (a^p)%c:
  242. //
  243. template <class Backend>
  244. BOOST_MP_CXX14_CONSTEXPR void eval_powm(Backend& result, const Backend& a, const Backend& p, const Backend& c)
  245. {
  246. using default_ops::eval_bit_test;
  247. using default_ops::eval_get_sign;
  248. using default_ops::eval_modulus;
  249. using default_ops::eval_multiply;
  250. using default_ops::eval_right_shift;
  251. using double_type = typename double_precision_type<Backend>::type ;
  252. using ui_type = typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type;
  253. check_sign_of_backend(p, std::integral_constant<bool, std::numeric_limits<number<Backend> >::is_signed>());
  254. double_type x, y(a), b(p), t;
  255. x = ui_type(1u);
  256. while (eval_get_sign(b) > 0)
  257. {
  258. if (eval_bit_test(b, 0))
  259. {
  260. eval_multiply(t, x, y);
  261. eval_modulus(x, t, c);
  262. }
  263. eval_multiply(t, y, y);
  264. eval_modulus(y, t, c);
  265. eval_right_shift(b, ui_type(1));
  266. }
  267. Backend x2(x);
  268. eval_modulus(result, x2, c);
  269. }
  270. template <class Backend, class Integer>
  271. BOOST_MP_CXX14_CONSTEXPR void eval_powm(Backend& result, const Backend& a, const Backend& p, Integer c)
  272. {
  273. using double_type = typename double_precision_type<Backend>::type ;
  274. using ui_type = typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type;
  275. using i1_type = typename boost::multiprecision::detail::canonical<Integer, double_type>::type ;
  276. using i2_type = typename boost::multiprecision::detail::canonical<Integer, Backend>::type ;
  277. using default_ops::eval_bit_test;
  278. using default_ops::eval_get_sign;
  279. using default_ops::eval_modulus;
  280. using default_ops::eval_multiply;
  281. using default_ops::eval_right_shift;
  282. check_sign_of_backend(p, std::integral_constant<bool, std::numeric_limits<number<Backend> >::is_signed>());
  283. if (eval_get_sign(p) < 0)
  284. {
  285. BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
  286. }
  287. double_type x, y(a), b(p), t;
  288. x = ui_type(1u);
  289. while (eval_get_sign(b) > 0)
  290. {
  291. if (eval_bit_test(b, 0))
  292. {
  293. eval_multiply(t, x, y);
  294. eval_modulus(x, t, static_cast<i1_type>(c));
  295. }
  296. eval_multiply(t, y, y);
  297. eval_modulus(y, t, static_cast<i1_type>(c));
  298. eval_right_shift(b, ui_type(1));
  299. }
  300. Backend x2(x);
  301. eval_modulus(result, x2, static_cast<i2_type>(c));
  302. }
  303. template <class Backend, class Integer>
  304. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_unsigned<Integer>::value >::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
  305. {
  306. using double_type = typename double_precision_type<Backend>::type ;
  307. using ui_type = typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type;
  308. using default_ops::eval_bit_test;
  309. using default_ops::eval_get_sign;
  310. using default_ops::eval_modulus;
  311. using default_ops::eval_multiply;
  312. using default_ops::eval_right_shift;
  313. double_type x, y(a), t;
  314. x = ui_type(1u);
  315. while (b > 0)
  316. {
  317. if (b & 1)
  318. {
  319. eval_multiply(t, x, y);
  320. eval_modulus(x, t, c);
  321. }
  322. eval_multiply(t, y, y);
  323. eval_modulus(y, t, c);
  324. b >>= 1;
  325. }
  326. Backend x2(x);
  327. eval_modulus(result, x2, c);
  328. }
  329. template <class Backend, class Integer>
  330. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_signed<Integer>::value && boost::multiprecision::detail::is_integral<Integer>::value>::type eval_powm(Backend& result, const Backend& a, Integer b, const Backend& c)
  331. {
  332. if (b < 0)
  333. {
  334. BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
  335. }
  336. eval_powm(result, a, static_cast<typename boost::multiprecision::detail::make_unsigned<Integer>::type>(b), c);
  337. }
  338. template <class Backend, class Integer1, class Integer2>
  339. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_unsigned<Integer1>::value >::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
  340. {
  341. using double_type = typename double_precision_type<Backend>::type ;
  342. using ui_type = typename boost::multiprecision::detail::canonical<unsigned char, double_type>::type;
  343. using i1_type = typename boost::multiprecision::detail::canonical<Integer1, double_type>::type ;
  344. using i2_type = typename boost::multiprecision::detail::canonical<Integer2, Backend>::type ;
  345. using default_ops::eval_bit_test;
  346. using default_ops::eval_get_sign;
  347. using default_ops::eval_modulus;
  348. using default_ops::eval_multiply;
  349. using default_ops::eval_right_shift;
  350. double_type x, y(a), t;
  351. x = ui_type(1u);
  352. while (b > 0)
  353. {
  354. if (b & 1)
  355. {
  356. eval_multiply(t, x, y);
  357. eval_modulus(x, t, static_cast<i1_type>(c));
  358. }
  359. eval_multiply(t, y, y);
  360. eval_modulus(y, t, static_cast<i1_type>(c));
  361. b >>= 1;
  362. }
  363. Backend x2(x);
  364. eval_modulus(result, x2, static_cast<i2_type>(c));
  365. }
  366. template <class Backend, class Integer1, class Integer2>
  367. BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_signed<Integer1>::value && boost::multiprecision::detail::is_integral<Integer1>::value>::type eval_powm(Backend& result, const Backend& a, Integer1 b, Integer2 c)
  368. {
  369. if (b < 0)
  370. {
  371. BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent."));
  372. }
  373. eval_powm(result, a, static_cast<typename boost::multiprecision::detail::make_unsigned<Integer1>::type>(b), c);
  374. }
  375. struct powm_func
  376. {
  377. template <class T, class U, class V>
  378. BOOST_MP_CXX14_CONSTEXPR void operator()(T& result, const T& b, const U& p, const V& m) const
  379. {
  380. eval_powm(result, b, p, m);
  381. }
  382. template <class R, class T, class U, class V>
  383. BOOST_MP_CXX14_CONSTEXPR void operator()(R& result, const T& b, const U& p, const V& m) const
  384. {
  385. T temp;
  386. eval_powm(temp, b, p, m);
  387. result = std::move(temp);
  388. }
  389. };
  390. } // namespace default_ops
  391. template <class T, class U, class V>
  392. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  393. (number_category<T>::value == number_kind_integer) &&
  394. (is_number<T>::value || is_number_expression<T>::value) &&
  395. (is_number<U>::value || is_number_expression<U>::value || boost::multiprecision::detail::is_integral<U>::value) &&
  396. (is_number<V>::value || is_number_expression<V>::value || boost::multiprecision::detail::is_integral<V>::value),
  397. typename std::conditional<
  398. is_no_et_number<T>::value,
  399. T,
  400. typename std::conditional<
  401. is_no_et_number<U>::value,
  402. U,
  403. typename std::conditional<
  404. is_no_et_number<V>::value,
  405. V,
  406. detail::expression<detail::function, default_ops::powm_func, T, U, V> >::type>::type>::type>::type
  407. powm(const T& b, const U& p, const V& mod)
  408. {
  409. return detail::expression<detail::function, default_ops::powm_func, T, U, V>(
  410. default_ops::powm_func(), b, p, mod);
  411. }
  412. }} // namespace boost::multiprecision
  413. #endif