integer_ops.hpp 18 KB

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