integer_ops.hpp 17 KB

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