generic_interconvert.hpp 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464
  1. ///////////////////////////////////////////////////////////////////////////////
  2. // Copyright 2011 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_0.txt)
  5. #ifndef BOOST_MP_GENERIC_INTERCONVERT_HPP
  6. #define BOOST_MP_GENERIC_INTERCONVERT_HPP
  7. #include <boost/multiprecision/detail/default_ops.hpp>
  8. #ifdef BOOST_MSVC
  9. #pragma warning(push)
  10. #pragma warning(disable:4127 6326)
  11. #endif
  12. namespace boost{ namespace multiprecision{ namespace detail{
  13. template <class To, class From>
  14. inline To do_cast(const From & from)
  15. {
  16. return static_cast<To>(from);
  17. }
  18. template <class To, class B, ::boost::multiprecision::expression_template_option et>
  19. inline To do_cast(const number<B, et>& from)
  20. {
  21. return from.template convert_to<To>();
  22. }
  23. template <class To, class From>
  24. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_floating_point>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/)
  25. {
  26. using default_ops::eval_get_sign;
  27. using default_ops::eval_bitwise_and;
  28. using default_ops::eval_convert_to;
  29. using default_ops::eval_right_shift;
  30. using default_ops::eval_ldexp;
  31. using default_ops::eval_add;
  32. // smallest unsigned type handled natively by "From" is likely to be it's limb_type:
  33. typedef typename canonical<unsigned char, From>::type limb_type;
  34. // get the corresponding type that we can assign to "To":
  35. typedef typename canonical<limb_type, To>::type to_type;
  36. From t(from);
  37. bool is_neg = eval_get_sign(t) < 0;
  38. if(is_neg)
  39. t.negate();
  40. // Pick off the first limb:
  41. limb_type limb;
  42. limb_type mask = ~static_cast<limb_type>(0);
  43. From fl;
  44. eval_bitwise_and(fl, t, mask);
  45. eval_convert_to(&limb, fl);
  46. to = static_cast<to_type>(limb);
  47. eval_right_shift(t, std::numeric_limits<limb_type>::digits);
  48. //
  49. // Then keep picking off more limbs until "t" is zero:
  50. //
  51. To l;
  52. unsigned shift = std::numeric_limits<limb_type>::digits;
  53. while(!eval_is_zero(t))
  54. {
  55. eval_bitwise_and(fl, t, mask);
  56. eval_convert_to(&limb, fl);
  57. l = static_cast<to_type>(limb);
  58. eval_right_shift(t, std::numeric_limits<limb_type>::digits);
  59. eval_ldexp(l, l, shift);
  60. eval_add(to, l);
  61. shift += std::numeric_limits<limb_type>::digits;
  62. }
  63. //
  64. // Finish off by setting the sign:
  65. //
  66. if(is_neg)
  67. to.negate();
  68. }
  69. template <class To, class From>
  70. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/)
  71. {
  72. using default_ops::eval_get_sign;
  73. using default_ops::eval_bitwise_and;
  74. using default_ops::eval_convert_to;
  75. using default_ops::eval_right_shift;
  76. using default_ops::eval_left_shift;
  77. using default_ops::eval_bitwise_or;
  78. using default_ops::eval_is_zero;
  79. // smallest unsigned type handled natively by "From" is likely to be it's limb_type:
  80. typedef typename canonical<unsigned char, From>::type limb_type;
  81. // get the corresponding type that we can assign to "To":
  82. typedef typename canonical<limb_type, To>::type to_type;
  83. From t(from);
  84. bool is_neg = eval_get_sign(t) < 0;
  85. if(is_neg)
  86. t.negate();
  87. // Pick off the first limb:
  88. limb_type limb;
  89. limb_type mask = static_cast<limb_type>(~static_cast<limb_type>(0));
  90. From fl;
  91. eval_bitwise_and(fl, t, mask);
  92. eval_convert_to(&limb, fl);
  93. to = static_cast<to_type>(limb);
  94. eval_right_shift(t, std::numeric_limits<limb_type>::digits);
  95. //
  96. // Then keep picking off more limbs until "t" is zero:
  97. //
  98. To l;
  99. unsigned shift = std::numeric_limits<limb_type>::digits;
  100. while(!eval_is_zero(t))
  101. {
  102. eval_bitwise_and(fl, t, mask);
  103. eval_convert_to(&limb, fl);
  104. l = static_cast<to_type>(limb);
  105. eval_right_shift(t, std::numeric_limits<limb_type>::digits);
  106. eval_left_shift(l, shift);
  107. eval_bitwise_or(to, l);
  108. shift += std::numeric_limits<limb_type>::digits;
  109. }
  110. //
  111. // Finish off by setting the sign:
  112. //
  113. if(is_neg)
  114. to.negate();
  115. }
  116. template <class To, class From>
  117. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_floating_point>& /*to_type*/, const mpl::int_<number_kind_floating_point>& /*from_type*/)
  118. {
  119. #ifdef BOOST_MSVC
  120. #pragma warning(push)
  121. #pragma warning(disable:4127)
  122. #endif
  123. //
  124. // The code here only works when the radix of "From" is 2, we could try shifting by other
  125. // radixes but it would complicate things.... use a string conversion when the radix is other
  126. // than 2:
  127. //
  128. if(std::numeric_limits<number<From> >::radix != 2)
  129. {
  130. to = from.str(0, std::ios_base::fmtflags()).c_str();
  131. return;
  132. }
  133. typedef typename canonical<unsigned char, To>::type ui_type;
  134. using default_ops::eval_fpclassify;
  135. using default_ops::eval_add;
  136. using default_ops::eval_subtract;
  137. using default_ops::eval_convert_to;
  138. //
  139. // First classify the input, then handle the special cases:
  140. //
  141. int c = eval_fpclassify(from);
  142. if(c == (int)FP_ZERO)
  143. {
  144. to = ui_type(0);
  145. return;
  146. }
  147. else if(c == (int)FP_NAN)
  148. {
  149. to = static_cast<const char*>("nan");
  150. return;
  151. }
  152. else if(c == (int)FP_INFINITE)
  153. {
  154. to = static_cast<const char*>("inf");
  155. if(eval_get_sign(from) < 0)
  156. to.negate();
  157. return;
  158. }
  159. typename From::exponent_type e;
  160. From f, term;
  161. to = ui_type(0);
  162. eval_frexp(f, from, &e);
  163. static const int shift = std::numeric_limits<boost::intmax_t>::digits - 1;
  164. while(!eval_is_zero(f))
  165. {
  166. // extract int sized bits from f:
  167. eval_ldexp(f, f, shift);
  168. eval_floor(term, f);
  169. e -= shift;
  170. eval_ldexp(to, to, shift);
  171. typename boost::multiprecision::detail::canonical<boost::intmax_t, To>::type ll;
  172. eval_convert_to(&ll, term);
  173. eval_add(to, ll);
  174. eval_subtract(f, term);
  175. }
  176. typedef typename To::exponent_type to_exponent;
  177. if((e > (std::numeric_limits<to_exponent>::max)()) || (e < (std::numeric_limits<to_exponent>::min)()))
  178. {
  179. to = static_cast<const char*>("inf");
  180. if(eval_get_sign(from) < 0)
  181. to.negate();
  182. return;
  183. }
  184. eval_ldexp(to, to, static_cast<to_exponent>(e));
  185. #ifdef BOOST_MSVC
  186. #pragma warning(pop)
  187. #endif
  188. }
  189. template <class To, class From>
  190. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/)
  191. {
  192. typedef typename component_type<number<To> >::type to_component_type;
  193. number<From> t(from);
  194. to_component_type n(numerator(t)), d(denominator(t));
  195. using default_ops::assign_components;
  196. assign_components(to, n.backend(), d.backend());
  197. }
  198. template <class To, class From>
  199. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/)
  200. {
  201. typedef typename component_type<number<To> >::type to_component_type;
  202. number<From> t(from);
  203. to_component_type n(t), d(1);
  204. using default_ops::assign_components;
  205. assign_components(to, n.backend(), d.backend());
  206. }
  207. template <class R, class LargeInteger>
  208. R safe_convert_to_float(const LargeInteger& i)
  209. {
  210. using std::ldexp;
  211. if(!i)
  212. return R(0);
  213. if(std::numeric_limits<R>::is_specialized && std::numeric_limits<R>::max_exponent)
  214. {
  215. LargeInteger val(i);
  216. if(val.sign() < 0)
  217. val = -val;
  218. unsigned mb = msb(val);
  219. if(mb >= std::numeric_limits<R>::max_exponent)
  220. {
  221. int scale_factor = (int)mb + 1 - std::numeric_limits<R>::max_exponent;
  222. BOOST_ASSERT(scale_factor >= 1);
  223. val >>= scale_factor;
  224. R result = val.template convert_to<R>();
  225. if(std::numeric_limits<R>::digits == 0 || std::numeric_limits<R>::digits >= std::numeric_limits<R>::max_exponent)
  226. {
  227. //
  228. // Calculate and add on the remainder, only if there are more
  229. // digits in the mantissa that the size of the exponent, in
  230. // other words if we are dropping digits in the conversion
  231. // otherwise:
  232. //
  233. LargeInteger remainder(i);
  234. remainder &= (LargeInteger(1) << scale_factor) - 1;
  235. result += ldexp(safe_convert_to_float<R>(remainder), -scale_factor);
  236. }
  237. return i.sign() < 0 ? static_cast<R>(-result) : result;
  238. }
  239. }
  240. return i.template convert_to<R>();
  241. }
  242. template <class To, class Integer>
  243. inline typename disable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  244. generic_convert_rational_to_float_imp(To& result, const Integer& n, const Integer& d, const mpl::true_&)
  245. {
  246. //
  247. // If we get here, then there's something about one type or the other
  248. // that prevents an exactly rounded result from being calculated
  249. // (or at least it's not clear how to implement such a thing).
  250. //
  251. using default_ops::eval_divide;
  252. number<To> fn(safe_convert_to_float<number<To> >(n)), fd(safe_convert_to_float<number<To> >(d));
  253. eval_divide(result, fn.backend(), fd.backend());
  254. }
  255. template <class To, class Integer>
  256. inline typename enable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  257. generic_convert_rational_to_float_imp(To& result, const Integer& n, const Integer& d, const mpl::true_&)
  258. {
  259. //
  260. // If we get here, then there's something about one type or the other
  261. // that prevents an exactly rounded result from being calculated
  262. // (or at least it's not clear how to implement such a thing).
  263. //
  264. To fd(safe_convert_to_float<To>(d));
  265. result = safe_convert_to_float<To>(n);
  266. result /= fd;
  267. }
  268. template <class To, class Integer>
  269. typename enable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  270. generic_convert_rational_to_float_imp(To& result, Integer& num, Integer& denom, const mpl::false_&)
  271. {
  272. //
  273. // If we get here, then the precision of type To is known, and the integer type is unbounded
  274. // so we can use integer division plus manipulation of the remainder to get an exactly
  275. // rounded result.
  276. //
  277. if(num == 0)
  278. {
  279. result = 0;
  280. return;
  281. }
  282. bool s = false;
  283. if(num < 0)
  284. {
  285. s = true;
  286. num = -num;
  287. }
  288. int denom_bits = msb(denom);
  289. int shift = std::numeric_limits<To>::digits + denom_bits - msb(num) + 1;
  290. if(shift > 0)
  291. num <<= shift;
  292. else if(shift < 0)
  293. denom <<= boost::multiprecision::detail::unsigned_abs(shift);
  294. Integer q, r;
  295. divide_qr(num, denom, q, r);
  296. int q_bits = msb(q);
  297. if(q_bits == std::numeric_limits<To>::digits)
  298. {
  299. //
  300. // Round up if 2 * r > denom:
  301. //
  302. r <<= 1;
  303. int c = r.compare(denom);
  304. if(c > 0)
  305. ++q;
  306. else if((c == 0) && (q & 1u))
  307. {
  308. ++q;
  309. }
  310. }
  311. else
  312. {
  313. BOOST_ASSERT(q_bits == 1 + std::numeric_limits<To>::digits);
  314. //
  315. // We basically already have the rounding info:
  316. //
  317. if(q & 1u)
  318. {
  319. if(r || (q & 2u))
  320. ++q;
  321. }
  322. }
  323. using std::ldexp;
  324. result = do_cast<To>(q);
  325. result = ldexp(result, -shift);
  326. if(s)
  327. result = -result;
  328. }
  329. template <class To, class Integer>
  330. inline typename disable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  331. generic_convert_rational_to_float_imp(To& result, Integer& num, Integer& denom, const mpl::false_& tag)
  332. {
  333. number<To> t;
  334. generic_convert_rational_to_float_imp(t, num, denom, tag);
  335. result = t.backend();
  336. }
  337. template <class To, class From>
  338. inline void generic_convert_rational_to_float(To& result, const From& f)
  339. {
  340. //
  341. // Type From is always a Backend to number<>, or an
  342. // instance of number<>, but we allow
  343. // To to be either a Backend type, or a real number type,
  344. // that way we can call this from generic conversions, and
  345. // from specific conversions to built in types.
  346. //
  347. typedef typename mpl::if_c<is_number<From>::value, From, number<From> >::type actual_from_type;
  348. typedef typename mpl::if_c<is_number<To>::value || is_floating_point<To>::value, To, number<To> >::type actual_to_type;
  349. typedef typename component_type<actual_from_type>::type integer_type;
  350. typedef mpl::bool_<!std::numeric_limits<integer_type>::is_specialized
  351. || std::numeric_limits<integer_type>::is_bounded
  352. || !std::numeric_limits<actual_to_type>::is_specialized
  353. || !std::numeric_limits<actual_to_type>::is_bounded
  354. || (std::numeric_limits<actual_to_type>::radix != 2)> dispatch_tag;
  355. integer_type n(numerator(static_cast<actual_from_type>(f))), d(denominator(static_cast<actual_from_type>(f)));
  356. generic_convert_rational_to_float_imp(result, n, d, dispatch_tag());
  357. }
  358. template <class To, class From>
  359. inline void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_floating_point>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/)
  360. {
  361. generic_convert_rational_to_float(to, from);
  362. }
  363. template <class To, class From>
  364. void generic_interconvert_float2rational(To& to, const From& from, const mpl::int_<2>& /*radix*/)
  365. {
  366. typedef typename mpl::front<typename To::unsigned_types>::type ui_type;
  367. static const int shift = std::numeric_limits<boost::long_long_type>::digits;
  368. typename From::exponent_type e;
  369. typename component_type<number<To> >::type num, denom;
  370. number<From> val(from);
  371. val = frexp(val, &e);
  372. while(val)
  373. {
  374. val = ldexp(val, shift);
  375. e -= shift;
  376. boost::long_long_type ll = boost::math::lltrunc(val);
  377. val -= ll;
  378. num <<= shift;
  379. num += ll;
  380. }
  381. denom = ui_type(1u);
  382. if(e < 0)
  383. denom <<= -e;
  384. else if(e > 0)
  385. num <<= e;
  386. assign_components(to, num.backend(), denom.backend());
  387. }
  388. template <class To, class From, int Radix>
  389. void generic_interconvert_float2rational(To& to, const From& from, const mpl::int_<Radix>& /*radix*/)
  390. {
  391. //
  392. // This is almost the same as the binary case above, but we have to use
  393. // scalbn and ilogb rather than ldexp and frexp, we also only extract
  394. // one Radix digit at a time which is terribly inefficient!
  395. //
  396. typedef typename mpl::front<typename To::unsigned_types>::type ui_type;
  397. typename From::exponent_type e;
  398. typename component_type<To>::type num, denom;
  399. number<From> val(from);
  400. e = ilogb(val);
  401. val = scalbn(val, -e);
  402. while(val)
  403. {
  404. boost::long_long_type ll = boost::math::lltrunc(val);
  405. val -= ll;
  406. val = scalbn(val, 1);
  407. num *= Radix;
  408. num += ll;
  409. --e;
  410. }
  411. ++e;
  412. denom = ui_type(Radix);
  413. denom = pow(denom, abs(e));
  414. if(e > 0)
  415. {
  416. num *= denom;
  417. denom = 1;
  418. }
  419. assign_components(to, num, denom);
  420. }
  421. template <class To, class From>
  422. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_floating_point>& /*from_type*/)
  423. {
  424. generic_interconvert_float2rational(to, from, mpl::int_<std::numeric_limits<number<From> >::radix>());
  425. }
  426. }}} // namespaces
  427. #ifdef BOOST_MSVC
  428. #pragma warning(pop)
  429. #endif
  430. #endif // BOOST_MP_GENERIC_INTERCONVERT_HPP