integer_ops.hpp 15 KB

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