generic_interconvert.hpp 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599
  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_add;
  27. using default_ops::eval_bitwise_and;
  28. using default_ops::eval_convert_to;
  29. using default_ops::eval_get_sign;
  30. using default_ops::eval_is_zero;
  31. using default_ops::eval_ldexp;
  32. using default_ops::eval_right_shift;
  33. // smallest unsigned type handled natively by "From" is likely to be it's limb_type:
  34. typedef typename canonical<unsigned char, From>::type l_limb_type;
  35. // get the corresponding type that we can assign to "To":
  36. typedef typename canonical<l_limb_type, To>::type to_type;
  37. From t(from);
  38. bool is_neg = eval_get_sign(t) < 0;
  39. if (is_neg)
  40. t.negate();
  41. // Pick off the first limb:
  42. l_limb_type limb;
  43. l_limb_type mask = static_cast<l_limb_type>(~static_cast<l_limb_type>(0));
  44. From fl;
  45. eval_bitwise_and(fl, t, mask);
  46. eval_convert_to(&limb, fl);
  47. to = static_cast<to_type>(limb);
  48. eval_right_shift(t, std::numeric_limits<l_limb_type>::digits);
  49. //
  50. // Then keep picking off more limbs until "t" is zero:
  51. //
  52. To l;
  53. unsigned shift = std::numeric_limits<l_limb_type>::digits;
  54. while (!eval_is_zero(t))
  55. {
  56. eval_bitwise_and(fl, t, mask);
  57. eval_convert_to(&limb, fl);
  58. l = static_cast<to_type>(limb);
  59. eval_right_shift(t, std::numeric_limits<l_limb_type>::digits);
  60. eval_ldexp(l, l, shift);
  61. eval_add(to, l);
  62. shift += std::numeric_limits<l_limb_type>::digits;
  63. }
  64. //
  65. // Finish off by setting the sign:
  66. //
  67. if (is_neg)
  68. to.negate();
  69. }
  70. template <class To, class From>
  71. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/)
  72. {
  73. using default_ops::eval_bitwise_and;
  74. using default_ops::eval_bitwise_or;
  75. using default_ops::eval_convert_to;
  76. using default_ops::eval_get_sign;
  77. using default_ops::eval_is_zero;
  78. using default_ops::eval_left_shift;
  79. using default_ops::eval_right_shift;
  80. // smallest unsigned type handled natively by "From" is likely to be it's limb_type:
  81. typedef typename canonical<unsigned char, From>::type limb_type;
  82. // get the corresponding type that we can assign to "To":
  83. typedef typename canonical<limb_type, To>::type to_type;
  84. From t(from);
  85. bool is_neg = eval_get_sign(t) < 0;
  86. if (is_neg)
  87. t.negate();
  88. // Pick off the first limb:
  89. limb_type limb;
  90. limb_type mask = static_cast<limb_type>(~static_cast<limb_type>(0));
  91. From fl;
  92. eval_bitwise_and(fl, t, mask);
  93. eval_convert_to(&limb, fl);
  94. to = static_cast<to_type>(limb);
  95. eval_right_shift(t, std::numeric_limits<limb_type>::digits);
  96. //
  97. // Then keep picking off more limbs until "t" is zero:
  98. //
  99. To l;
  100. unsigned shift = std::numeric_limits<limb_type>::digits;
  101. while (!eval_is_zero(t))
  102. {
  103. eval_bitwise_and(fl, t, mask);
  104. eval_convert_to(&limb, fl);
  105. l = static_cast<to_type>(limb);
  106. eval_right_shift(t, std::numeric_limits<limb_type>::digits);
  107. eval_left_shift(l, shift);
  108. eval_bitwise_or(to, l);
  109. shift += std::numeric_limits<limb_type>::digits;
  110. }
  111. //
  112. // Finish off by setting the sign:
  113. //
  114. if (is_neg)
  115. to.negate();
  116. }
  117. template <class To, class From>
  118. 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*/)
  119. {
  120. #ifdef BOOST_MSVC
  121. #pragma warning(push)
  122. #pragma warning(disable : 4127)
  123. #endif
  124. //
  125. // The code here only works when the radix of "From" is 2, we could try shifting by other
  126. // radixes but it would complicate things.... use a string conversion when the radix is other
  127. // than 2:
  128. //
  129. if (std::numeric_limits<number<From> >::radix != 2)
  130. {
  131. to = from.str(0, std::ios_base::fmtflags()).c_str();
  132. return;
  133. }
  134. typedef typename canonical<unsigned char, To>::type ui_type;
  135. using default_ops::eval_add;
  136. using default_ops::eval_convert_to;
  137. using default_ops::eval_fpclassify;
  138. using default_ops::eval_get_sign;
  139. using default_ops::eval_is_zero;
  140. using default_ops::eval_subtract;
  141. //
  142. // First classify the input, then handle the special cases:
  143. //
  144. int c = eval_fpclassify(from);
  145. if (c == (int)FP_ZERO)
  146. {
  147. to = ui_type(0);
  148. return;
  149. }
  150. else if (c == (int)FP_NAN)
  151. {
  152. to = static_cast<const char*>("nan");
  153. return;
  154. }
  155. else if (c == (int)FP_INFINITE)
  156. {
  157. to = static_cast<const char*>("inf");
  158. if (eval_get_sign(from) < 0)
  159. to.negate();
  160. return;
  161. }
  162. typename From::exponent_type e;
  163. From f, term;
  164. to = ui_type(0);
  165. eval_frexp(f, from, &e);
  166. static const int shift = std::numeric_limits<boost::intmax_t>::digits - 1;
  167. while (!eval_is_zero(f))
  168. {
  169. // extract int sized bits from f:
  170. eval_ldexp(f, f, shift);
  171. eval_floor(term, f);
  172. e -= shift;
  173. eval_ldexp(to, to, shift);
  174. typename boost::multiprecision::detail::canonical<boost::intmax_t, To>::type ll;
  175. eval_convert_to(&ll, term);
  176. eval_add(to, ll);
  177. eval_subtract(f, term);
  178. }
  179. typedef typename To::exponent_type to_exponent;
  180. if (e > (std::numeric_limits<to_exponent>::max)())
  181. {
  182. to = static_cast<const char*>("inf");
  183. if (eval_get_sign(from) < 0)
  184. to.negate();
  185. return;
  186. }
  187. if (e < (std::numeric_limits<to_exponent>::min)())
  188. {
  189. to = ui_type(0);
  190. if (eval_get_sign(from) < 0)
  191. to.negate();
  192. return;
  193. }
  194. eval_ldexp(to, to, static_cast<to_exponent>(e));
  195. #ifdef BOOST_MSVC
  196. #pragma warning(pop)
  197. #endif
  198. }
  199. template <class To, class From>
  200. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/)
  201. {
  202. typedef typename component_type<number<To> >::type to_component_type;
  203. number<From> t(from);
  204. to_component_type n(numerator(t)), d(denominator(t));
  205. using default_ops::assign_components;
  206. assign_components(to, n.backend(), d.backend());
  207. }
  208. template <class To, class From>
  209. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/)
  210. {
  211. typedef typename component_type<number<To> >::type to_component_type;
  212. number<From> t(from);
  213. to_component_type n(t), d(1);
  214. using default_ops::assign_components;
  215. assign_components(to, n.backend(), d.backend());
  216. }
  217. template <class R, class LargeInteger>
  218. R safe_convert_to_float(const LargeInteger& i)
  219. {
  220. using std::ldexp;
  221. if (!i)
  222. return R(0);
  223. if (std::numeric_limits<R>::is_specialized && std::numeric_limits<R>::max_exponent)
  224. {
  225. LargeInteger val(i);
  226. if (val.sign() < 0)
  227. val = -val;
  228. unsigned mb = msb(val);
  229. if (mb >= std::numeric_limits<R>::max_exponent)
  230. {
  231. int scale_factor = (int)mb + 1 - std::numeric_limits<R>::max_exponent;
  232. BOOST_ASSERT(scale_factor >= 1);
  233. val >>= scale_factor;
  234. R result = val.template convert_to<R>();
  235. if (std::numeric_limits<R>::digits == 0 || std::numeric_limits<R>::digits >= std::numeric_limits<R>::max_exponent)
  236. {
  237. //
  238. // Calculate and add on the remainder, only if there are more
  239. // digits in the mantissa that the size of the exponent, in
  240. // other words if we are dropping digits in the conversion
  241. // otherwise:
  242. //
  243. LargeInteger remainder(i);
  244. remainder &= (LargeInteger(1) << scale_factor) - 1;
  245. result += ldexp(safe_convert_to_float<R>(remainder), -scale_factor);
  246. }
  247. return i.sign() < 0 ? static_cast<R>(-result) : result;
  248. }
  249. }
  250. return i.template convert_to<R>();
  251. }
  252. template <class To, class Integer>
  253. inline typename disable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  254. generic_convert_rational_to_float_imp(To& result, const Integer& n, const Integer& d, const mpl::true_&)
  255. {
  256. //
  257. // If we get here, then there's something about one type or the other
  258. // that prevents an exactly rounded result from being calculated
  259. // (or at least it's not clear how to implement such a thing).
  260. //
  261. using default_ops::eval_divide;
  262. number<To> fn(safe_convert_to_float<number<To> >(n)), fd(safe_convert_to_float<number<To> >(d));
  263. eval_divide(result, fn.backend(), fd.backend());
  264. }
  265. template <class To, class Integer>
  266. inline typename enable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  267. generic_convert_rational_to_float_imp(To& result, const Integer& n, const Integer& d, const mpl::true_&)
  268. {
  269. //
  270. // If we get here, then there's something about one type or the other
  271. // that prevents an exactly rounded result from being calculated
  272. // (or at least it's not clear how to implement such a thing).
  273. //
  274. To fd(safe_convert_to_float<To>(d));
  275. result = safe_convert_to_float<To>(n);
  276. result /= fd;
  277. }
  278. template <class To, class Integer>
  279. typename enable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  280. generic_convert_rational_to_float_imp(To& result, Integer& num, Integer& denom, const mpl::false_&)
  281. {
  282. //
  283. // If we get here, then the precision of type To is known, and the integer type is unbounded
  284. // so we can use integer division plus manipulation of the remainder to get an exactly
  285. // rounded result.
  286. //
  287. if (num == 0)
  288. {
  289. result = 0;
  290. return;
  291. }
  292. bool s = false;
  293. if (num < 0)
  294. {
  295. s = true;
  296. num = -num;
  297. }
  298. int denom_bits = msb(denom);
  299. int shift = std::numeric_limits<To>::digits + denom_bits - msb(num);
  300. if (shift > 0)
  301. num <<= shift;
  302. else if (shift < 0)
  303. denom <<= boost::multiprecision::detail::unsigned_abs(shift);
  304. Integer q, r;
  305. divide_qr(num, denom, q, r);
  306. int q_bits = msb(q);
  307. if (q_bits == std::numeric_limits<To>::digits - 1)
  308. {
  309. //
  310. // Round up if 2 * r > denom:
  311. //
  312. r <<= 1;
  313. int c = r.compare(denom);
  314. if (c > 0)
  315. ++q;
  316. else if ((c == 0) && (q & 1u))
  317. {
  318. ++q;
  319. }
  320. }
  321. else
  322. {
  323. BOOST_ASSERT(q_bits == std::numeric_limits<To>::digits);
  324. //
  325. // We basically already have the rounding info:
  326. //
  327. if (q & 1u)
  328. {
  329. if (r || (q & 2u))
  330. ++q;
  331. }
  332. }
  333. using std::ldexp;
  334. result = do_cast<To>(q);
  335. result = ldexp(result, -shift);
  336. if (s)
  337. result = -result;
  338. }
  339. template <class To, class Integer>
  340. inline typename disable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  341. generic_convert_rational_to_float_imp(To& result, Integer& num, Integer& denom, const mpl::false_& tag)
  342. {
  343. number<To> t;
  344. generic_convert_rational_to_float_imp(t, num, denom, tag);
  345. result = t.backend();
  346. }
  347. template <class To, class From>
  348. inline void generic_convert_rational_to_float(To& result, const From& f)
  349. {
  350. //
  351. // Type From is always a Backend to number<>, or an
  352. // instance of number<>, but we allow
  353. // To to be either a Backend type, or a real number type,
  354. // that way we can call this from generic conversions, and
  355. // from specific conversions to built in types.
  356. //
  357. typedef typename mpl::if_c<is_number<From>::value, From, number<From> >::type actual_from_type;
  358. typedef typename mpl::if_c<is_number<To>::value || is_floating_point<To>::value, To, number<To> >::type actual_to_type;
  359. typedef typename component_type<actual_from_type>::type integer_type;
  360. typedef mpl::bool_<!std::numeric_limits<integer_type>::is_specialized || std::numeric_limits<integer_type>::is_bounded || !std::numeric_limits<actual_to_type>::is_specialized || !std::numeric_limits<actual_to_type>::is_bounded || (std::numeric_limits<actual_to_type>::radix != 2)> dispatch_tag;
  361. integer_type n(numerator(static_cast<actual_from_type>(f))), d(denominator(static_cast<actual_from_type>(f)));
  362. generic_convert_rational_to_float_imp(result, n, d, dispatch_tag());
  363. }
  364. template <class To, class From>
  365. 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*/)
  366. {
  367. generic_convert_rational_to_float(to, from);
  368. }
  369. template <class To, class From>
  370. void generic_interconvert_float2rational(To& to, const From& from, const mpl::int_<2>& /*radix*/)
  371. {
  372. typedef typename mpl::front<typename To::unsigned_types>::type ui_type;
  373. static const int shift = std::numeric_limits<boost::long_long_type>::digits;
  374. typename From::exponent_type e;
  375. typename component_type<number<To> >::type num, denom;
  376. number<From> val(from);
  377. val = frexp(val, &e);
  378. while (val)
  379. {
  380. val = ldexp(val, shift);
  381. e -= shift;
  382. boost::long_long_type ll = boost::math::lltrunc(val);
  383. val -= ll;
  384. num <<= shift;
  385. num += ll;
  386. }
  387. denom = ui_type(1u);
  388. if (e < 0)
  389. denom <<= -e;
  390. else if (e > 0)
  391. num <<= e;
  392. assign_components(to, num.backend(), denom.backend());
  393. }
  394. template <class To, class From, int Radix>
  395. void generic_interconvert_float2rational(To& to, const From& from, const mpl::int_<Radix>& /*radix*/)
  396. {
  397. //
  398. // This is almost the same as the binary case above, but we have to use
  399. // scalbn and ilogb rather than ldexp and frexp, we also only extract
  400. // one Radix digit at a time which is terribly inefficient!
  401. //
  402. typedef typename mpl::front<typename To::unsigned_types>::type ui_type;
  403. typename From::exponent_type e;
  404. typename component_type<number<To> >::type num, denom;
  405. number<From> val(from);
  406. if (!val)
  407. {
  408. to = ui_type(0u);
  409. return;
  410. }
  411. e = ilogb(val);
  412. val = scalbn(val, -e);
  413. while (val)
  414. {
  415. boost::long_long_type ll = boost::math::lltrunc(val);
  416. val -= ll;
  417. val = scalbn(val, 1);
  418. num *= Radix;
  419. num += ll;
  420. --e;
  421. }
  422. ++e;
  423. denom = ui_type(Radix);
  424. denom = pow(denom, abs(e));
  425. if (e > 0)
  426. {
  427. num *= denom;
  428. denom = 1;
  429. }
  430. assign_components(to, num.backend(), denom.backend());
  431. }
  432. template <class To, class From>
  433. 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*/)
  434. {
  435. generic_interconvert_float2rational(to, from, mpl::int_<std::numeric_limits<number<From> >::radix>());
  436. }
  437. template <class To, class From>
  438. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/)
  439. {
  440. number<From> t(from);
  441. number<To> result(numerator(t) / denominator(t));
  442. to = result.backend();
  443. }
  444. template <class To, class From>
  445. void generic_interconvert_float2int(To& to, const From& from, const mpl::int_<2>& /*radix*/)
  446. {
  447. typedef typename From::exponent_type exponent_type;
  448. static const exponent_type shift = std::numeric_limits<boost::long_long_type>::digits;
  449. exponent_type e;
  450. number<To> num(0u);
  451. number<From> val(from);
  452. val = frexp(val, &e);
  453. bool neg = false;
  454. if (val.sign() < 0)
  455. {
  456. val.backend().negate();
  457. neg = true;
  458. }
  459. while (e > 0)
  460. {
  461. exponent_type s = (std::min)(e, shift);
  462. val = ldexp(val, s);
  463. e -= s;
  464. boost::long_long_type ll = boost::math::lltrunc(val);
  465. val -= ll;
  466. num <<= s;
  467. num += ll;
  468. }
  469. to = num.backend();
  470. if (neg)
  471. to.negate();
  472. }
  473. template <class To, class From, int Radix>
  474. void generic_interconvert_float2int(To& to, const From& from, const mpl::int_<Radix>& /*radix*/)
  475. {
  476. //
  477. // This is almost the same as the binary case above, but we have to use
  478. // scalbn and ilogb rather than ldexp and frexp, we also only extract
  479. // one Radix digit at a time which is terribly inefficient!
  480. //
  481. typename From::exponent_type e;
  482. number<To> num(0u);
  483. number<From> val(from);
  484. e = ilogb(val);
  485. val = scalbn(val, -e);
  486. while (e >= 0)
  487. {
  488. boost::long_long_type ll = boost::math::lltrunc(val);
  489. val -= ll;
  490. val = scalbn(val, 1);
  491. num *= Radix;
  492. num += ll;
  493. --e;
  494. }
  495. to = num.backend();
  496. }
  497. template <class To, class From>
  498. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_floating_point>& /*from_type*/)
  499. {
  500. generic_interconvert_float2int(to, from, mpl::int_<std::numeric_limits<number<From> >::radix>());
  501. }
  502. template <class To, class From, class tag>
  503. void generic_interconvert_complex_to_scalar(To& to, const From& from, const mpl::true_&, const tag&)
  504. {
  505. // We just want the real part, and "to" is the correct type already:
  506. eval_real(to, from);
  507. To im;
  508. eval_imag(im, from);
  509. if (!eval_is_zero(im))
  510. BOOST_THROW_EXCEPTION(std::runtime_error("Could not convert imaginary number to scalar."));
  511. }
  512. template <class To, class From>
  513. void generic_interconvert_complex_to_scalar(To& to, const From& from, const mpl::false_&, const mpl::true_&)
  514. {
  515. typedef typename component_type<number<From> >::type component_number;
  516. typedef typename component_number::backend_type component_backend;
  517. //
  518. // Get the real part and copy-construct the result from it:
  519. //
  520. component_backend r;
  521. generic_interconvert_complex_to_scalar(r, from, mpl::true_(), mpl::true_());
  522. to = r;
  523. }
  524. template <class To, class From>
  525. void generic_interconvert_complex_to_scalar(To& to, const From& from, const mpl::false_&, const mpl::false_&)
  526. {
  527. typedef typename component_type<number<From> >::type component_number;
  528. typedef typename component_number::backend_type component_backend;
  529. //
  530. // Get the real part and use a generic_interconvert to type To:
  531. //
  532. component_backend r;
  533. generic_interconvert_complex_to_scalar(r, from, mpl::true_(), mpl::true_());
  534. generic_interconvert(to, r, mpl::int_<number_category<To>::value>(), mpl::int_<number_category<To>::value>());
  535. }
  536. template <class To, class From>
  537. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_floating_point>& /*to_type*/, const mpl::int_<number_kind_complex>& /*from_type*/)
  538. {
  539. typedef typename component_type<number<From> >::type component_number;
  540. typedef typename component_number::backend_type component_backend;
  541. generic_interconvert_complex_to_scalar(to, from, mpl::bool_<boost::is_same<component_backend, To>::value>(), mpl::bool_<boost::is_constructible<To, const component_backend&>::value>());
  542. }
  543. template <class To, class From>
  544. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_complex>& /*from_type*/)
  545. {
  546. typedef typename component_type<number<From> >::type component_number;
  547. typedef typename component_number::backend_type component_backend;
  548. generic_interconvert_complex_to_scalar(to, from, mpl::bool_<boost::is_same<component_backend, To>::value>(), mpl::bool_<boost::is_constructible<To, const component_backend&>::value>());
  549. }
  550. }
  551. }
  552. } // namespace boost::multiprecision::detail
  553. #ifdef BOOST_MSVC
  554. #pragma warning(pop)
  555. #endif
  556. #endif // BOOST_MP_GENERIC_INTERCONVERT_HPP