generic_interconvert.hpp 20 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_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. using default_ops::eval_is_zero;
  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_get_sign;
  74. using default_ops::eval_bitwise_and;
  75. using default_ops::eval_convert_to;
  76. using default_ops::eval_right_shift;
  77. using default_ops::eval_left_shift;
  78. using default_ops::eval_bitwise_or;
  79. using default_ops::eval_is_zero;
  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_fpclassify;
  136. using default_ops::eval_add;
  137. using default_ops::eval_subtract;
  138. using default_ops::eval_convert_to;
  139. using default_ops::eval_get_sign;
  140. using default_ops::eval_is_zero;
  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)()) || (e < (std::numeric_limits<to_exponent>::min)()))
  181. {
  182. to = static_cast<const char*>("inf");
  183. if(eval_get_sign(from) < 0)
  184. to.negate();
  185. return;
  186. }
  187. eval_ldexp(to, to, static_cast<to_exponent>(e));
  188. #ifdef BOOST_MSVC
  189. #pragma warning(pop)
  190. #endif
  191. }
  192. template <class To, class From>
  193. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/)
  194. {
  195. typedef typename component_type<number<To> >::type to_component_type;
  196. number<From> t(from);
  197. to_component_type n(numerator(t)), d(denominator(t));
  198. using default_ops::assign_components;
  199. assign_components(to, n.backend(), d.backend());
  200. }
  201. template <class To, class From>
  202. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/)
  203. {
  204. typedef typename component_type<number<To> >::type to_component_type;
  205. number<From> t(from);
  206. to_component_type n(t), d(1);
  207. using default_ops::assign_components;
  208. assign_components(to, n.backend(), d.backend());
  209. }
  210. template <class R, class LargeInteger>
  211. R safe_convert_to_float(const LargeInteger& i)
  212. {
  213. using std::ldexp;
  214. if(!i)
  215. return R(0);
  216. if(std::numeric_limits<R>::is_specialized && std::numeric_limits<R>::max_exponent)
  217. {
  218. LargeInteger val(i);
  219. if(val.sign() < 0)
  220. val = -val;
  221. unsigned mb = msb(val);
  222. if(mb >= std::numeric_limits<R>::max_exponent)
  223. {
  224. int scale_factor = (int)mb + 1 - std::numeric_limits<R>::max_exponent;
  225. BOOST_ASSERT(scale_factor >= 1);
  226. val >>= scale_factor;
  227. R result = val.template convert_to<R>();
  228. if(std::numeric_limits<R>::digits == 0 || std::numeric_limits<R>::digits >= std::numeric_limits<R>::max_exponent)
  229. {
  230. //
  231. // Calculate and add on the remainder, only if there are more
  232. // digits in the mantissa that the size of the exponent, in
  233. // other words if we are dropping digits in the conversion
  234. // otherwise:
  235. //
  236. LargeInteger remainder(i);
  237. remainder &= (LargeInteger(1) << scale_factor) - 1;
  238. result += ldexp(safe_convert_to_float<R>(remainder), -scale_factor);
  239. }
  240. return i.sign() < 0 ? static_cast<R>(-result) : result;
  241. }
  242. }
  243. return i.template convert_to<R>();
  244. }
  245. template <class To, class Integer>
  246. inline typename disable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  247. generic_convert_rational_to_float_imp(To& result, const Integer& n, const Integer& d, const mpl::true_&)
  248. {
  249. //
  250. // If we get here, then there's something about one type or the other
  251. // that prevents an exactly rounded result from being calculated
  252. // (or at least it's not clear how to implement such a thing).
  253. //
  254. using default_ops::eval_divide;
  255. number<To> fn(safe_convert_to_float<number<To> >(n)), fd(safe_convert_to_float<number<To> >(d));
  256. eval_divide(result, fn.backend(), fd.backend());
  257. }
  258. template <class To, class Integer>
  259. inline typename enable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  260. generic_convert_rational_to_float_imp(To& result, const Integer& n, const Integer& d, const mpl::true_&)
  261. {
  262. //
  263. // If we get here, then there's something about one type or the other
  264. // that prevents an exactly rounded result from being calculated
  265. // (or at least it's not clear how to implement such a thing).
  266. //
  267. To fd(safe_convert_to_float<To>(d));
  268. result = safe_convert_to_float<To>(n);
  269. result /= fd;
  270. }
  271. template <class To, class Integer>
  272. typename enable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  273. generic_convert_rational_to_float_imp(To& result, Integer& num, Integer& denom, const mpl::false_&)
  274. {
  275. //
  276. // If we get here, then the precision of type To is known, and the integer type is unbounded
  277. // so we can use integer division plus manipulation of the remainder to get an exactly
  278. // rounded result.
  279. //
  280. if(num == 0)
  281. {
  282. result = 0;
  283. return;
  284. }
  285. bool s = false;
  286. if(num < 0)
  287. {
  288. s = true;
  289. num = -num;
  290. }
  291. int denom_bits = msb(denom);
  292. int shift = std::numeric_limits<To>::digits + denom_bits - msb(num);
  293. if(shift > 0)
  294. num <<= shift;
  295. else if(shift < 0)
  296. denom <<= boost::multiprecision::detail::unsigned_abs(shift);
  297. Integer q, r;
  298. divide_qr(num, denom, q, r);
  299. int q_bits = msb(q);
  300. if(q_bits == std::numeric_limits<To>::digits - 1)
  301. {
  302. //
  303. // Round up if 2 * r > denom:
  304. //
  305. r <<= 1;
  306. int c = r.compare(denom);
  307. if(c > 0)
  308. ++q;
  309. else if((c == 0) && (q & 1u))
  310. {
  311. ++q;
  312. }
  313. }
  314. else
  315. {
  316. BOOST_ASSERT(q_bits == std::numeric_limits<To>::digits);
  317. //
  318. // We basically already have the rounding info:
  319. //
  320. if(q & 1u)
  321. {
  322. if(r || (q & 2u))
  323. ++q;
  324. }
  325. }
  326. using std::ldexp;
  327. result = do_cast<To>(q);
  328. result = ldexp(result, -shift);
  329. if(s)
  330. result = -result;
  331. }
  332. template <class To, class Integer>
  333. inline typename disable_if_c<is_number<To>::value || is_floating_point<To>::value>::type
  334. generic_convert_rational_to_float_imp(To& result, Integer& num, Integer& denom, const mpl::false_& tag)
  335. {
  336. number<To> t;
  337. generic_convert_rational_to_float_imp(t, num, denom, tag);
  338. result = t.backend();
  339. }
  340. template <class To, class From>
  341. inline void generic_convert_rational_to_float(To& result, const From& f)
  342. {
  343. //
  344. // Type From is always a Backend to number<>, or an
  345. // instance of number<>, but we allow
  346. // To to be either a Backend type, or a real number type,
  347. // that way we can call this from generic conversions, and
  348. // from specific conversions to built in types.
  349. //
  350. typedef typename mpl::if_c<is_number<From>::value, From, number<From> >::type actual_from_type;
  351. typedef typename mpl::if_c<is_number<To>::value || is_floating_point<To>::value, To, number<To> >::type actual_to_type;
  352. typedef typename component_type<actual_from_type>::type integer_type;
  353. typedef mpl::bool_<!std::numeric_limits<integer_type>::is_specialized
  354. || std::numeric_limits<integer_type>::is_bounded
  355. || !std::numeric_limits<actual_to_type>::is_specialized
  356. || !std::numeric_limits<actual_to_type>::is_bounded
  357. || (std::numeric_limits<actual_to_type>::radix != 2)> dispatch_tag;
  358. integer_type n(numerator(static_cast<actual_from_type>(f))), d(denominator(static_cast<actual_from_type>(f)));
  359. generic_convert_rational_to_float_imp(result, n, d, dispatch_tag());
  360. }
  361. template <class To, class From>
  362. 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*/)
  363. {
  364. generic_convert_rational_to_float(to, from);
  365. }
  366. template <class To, class From>
  367. void generic_interconvert_float2rational(To& to, const From& from, const mpl::int_<2>& /*radix*/)
  368. {
  369. typedef typename mpl::front<typename To::unsigned_types>::type ui_type;
  370. static const int shift = std::numeric_limits<boost::long_long_type>::digits;
  371. typename From::exponent_type e;
  372. typename component_type<number<To> >::type num, denom;
  373. number<From> val(from);
  374. val = frexp(val, &e);
  375. while(val)
  376. {
  377. val = ldexp(val, shift);
  378. e -= shift;
  379. boost::long_long_type ll = boost::math::lltrunc(val);
  380. val -= ll;
  381. num <<= shift;
  382. num += ll;
  383. }
  384. denom = ui_type(1u);
  385. if(e < 0)
  386. denom <<= -e;
  387. else if(e > 0)
  388. num <<= e;
  389. assign_components(to, num.backend(), denom.backend());
  390. }
  391. template <class To, class From, int Radix>
  392. void generic_interconvert_float2rational(To& to, const From& from, const mpl::int_<Radix>& /*radix*/)
  393. {
  394. //
  395. // This is almost the same as the binary case above, but we have to use
  396. // scalbn and ilogb rather than ldexp and frexp, we also only extract
  397. // one Radix digit at a time which is terribly inefficient!
  398. //
  399. typedef typename mpl::front<typename To::unsigned_types>::type ui_type;
  400. typename From::exponent_type e;
  401. typename component_type<number<To> >::type num, denom;
  402. number<From> val(from);
  403. if (!val)
  404. {
  405. to = ui_type(0u);
  406. return;
  407. }
  408. e = ilogb(val);
  409. val = scalbn(val, -e);
  410. while(val)
  411. {
  412. boost::long_long_type ll = boost::math::lltrunc(val);
  413. val -= ll;
  414. val = scalbn(val, 1);
  415. num *= Radix;
  416. num += ll;
  417. --e;
  418. }
  419. ++e;
  420. denom = ui_type(Radix);
  421. denom = pow(denom, abs(e));
  422. if(e > 0)
  423. {
  424. num *= denom;
  425. denom = 1;
  426. }
  427. assign_components(to, num.backend(), denom.backend());
  428. }
  429. template <class To, class From>
  430. 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*/)
  431. {
  432. generic_interconvert_float2rational(to, from, mpl::int_<std::numeric_limits<number<From> >::radix>());
  433. }
  434. template <class To, class From>
  435. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/)
  436. {
  437. number<From> t(from);
  438. number<To> result(numerator(t) / denominator(t));
  439. to = result.backend();
  440. }
  441. template <class To, class From>
  442. void generic_interconvert_float2int(To& to, const From& from, const mpl::int_<2>& /*radix*/)
  443. {
  444. typedef typename From::exponent_type exponent_type;
  445. static const exponent_type shift = std::numeric_limits<boost::long_long_type>::digits;
  446. exponent_type e;
  447. number<To> num(0u);
  448. number<From> val(from);
  449. val = frexp(val, &e);
  450. bool neg = false;
  451. if (val.sign() < 0)
  452. {
  453. val.backend().negate();
  454. neg = true;
  455. }
  456. while(e > 0)
  457. {
  458. exponent_type s = (std::min)(e, shift);
  459. val = ldexp(val, s);
  460. e -= s;
  461. boost::long_long_type ll = boost::math::lltrunc(val);
  462. val -= ll;
  463. num <<= s;
  464. num += ll;
  465. }
  466. to = num.backend();
  467. if (neg)
  468. to.negate();
  469. }
  470. template <class To, class From, int Radix>
  471. void generic_interconvert_float2int(To& to, const From& from, const mpl::int_<Radix>& /*radix*/)
  472. {
  473. //
  474. // This is almost the same as the binary case above, but we have to use
  475. // scalbn and ilogb rather than ldexp and frexp, we also only extract
  476. // one Radix digit at a time which is terribly inefficient!
  477. //
  478. typename From::exponent_type e;
  479. number<To> num(0u);
  480. number<From> val(from);
  481. e = ilogb(val);
  482. val = scalbn(val, -e);
  483. while(e >= 0)
  484. {
  485. boost::long_long_type ll = boost::math::lltrunc(val);
  486. val -= ll;
  487. val = scalbn(val, 1);
  488. num *= Radix;
  489. num += ll;
  490. --e;
  491. }
  492. to = num.backend();
  493. }
  494. template <class To, class From>
  495. 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*/)
  496. {
  497. generic_interconvert_float2int(to, from, mpl::int_<std::numeric_limits<number<From> >::radix>());
  498. }
  499. template <class To, class From, class tag>
  500. void generic_interconvert_complex_to_scalar(To& to, const From& from, const mpl::true_&, const tag&)
  501. {
  502. // We just want the real part, and "to" is the correct type already:
  503. eval_real(to, from);
  504. To im;
  505. eval_imag(im, from);
  506. if(!eval_is_zero(im))
  507. BOOST_THROW_EXCEPTION(std::runtime_error("Could not convert imaginary number to scalar."));
  508. }
  509. template <class To, class From>
  510. void generic_interconvert_complex_to_scalar(To& to, const From& from, const mpl::false_&, const mpl::true_&)
  511. {
  512. typedef typename component_type<number<From> >::type component_number;
  513. typedef typename component_number::backend_type component_backend;
  514. //
  515. // Get the real part and copy-construct the result from it:
  516. //
  517. component_backend r;
  518. generic_interconvert_complex_to_scalar(r, from, mpl::true_(), mpl::true_());
  519. to = r;
  520. }
  521. template <class To, class From>
  522. void generic_interconvert_complex_to_scalar(To& to, const From& from, const mpl::false_&, const mpl::false_&)
  523. {
  524. typedef typename component_type<number<From> >::type component_number;
  525. typedef typename component_number::backend_type component_backend;
  526. //
  527. // Get the real part and use a generic_interconvert to type To:
  528. //
  529. component_backend r;
  530. generic_interconvert_complex_to_scalar(r, from, mpl::true_(), mpl::true_());
  531. generic_interconvert(to, r, mpl::int_<number_category<To>::value>(), mpl::int_<number_category<To>::value>());
  532. }
  533. template <class To, class From>
  534. 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*/)
  535. {
  536. typedef typename component_type<number<From> >::type component_number;
  537. typedef typename component_number::backend_type component_backend;
  538. 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>());
  539. }
  540. template <class To, class From>
  541. void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_complex>& /*from_type*/)
  542. {
  543. typedef typename component_type<number<From> >::type component_number;
  544. typedef typename component_number::backend_type component_backend;
  545. 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>());
  546. }
  547. }
  548. }
  549. } // namespaces
  550. #ifdef BOOST_MSVC
  551. #pragma warning(pop)
  552. #endif
  553. #endif // BOOST_MP_GENERIC_INTERCONVERT_HPP