io.hpp 25 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697
  1. ///////////////////////////////////////////////////////////////
  2. // Copyright 2013 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_CPP_BIN_FLOAT_IO_HPP
  6. #define BOOST_MP_CPP_BIN_FLOAT_IO_HPP
  7. namespace boost { namespace multiprecision {
  8. namespace cpp_bf_io_detail {
  9. #ifdef BOOST_MSVC
  10. #pragma warning(push)
  11. #pragma warning(disable : 4127) // conditional expression is constant
  12. #endif
  13. //
  14. // Multiplies a by b and shifts the result so it fits inside max_bits bits,
  15. // returns by how much the result was shifted.
  16. //
  17. template <class I>
  18. inline I restricted_multiply(cpp_int& result, const cpp_int& a, const cpp_int& b, I max_bits, boost::int64_t& error)
  19. {
  20. result = a * b;
  21. I gb = msb(result);
  22. I rshift = 0;
  23. if (gb > max_bits)
  24. {
  25. rshift = gb - max_bits;
  26. I lb = lsb(result);
  27. int roundup = 0;
  28. // The error rate increases by the error of both a and b,
  29. // this may be overly pessimistic in many case as we're assuming
  30. // that a and b have the same level of uncertainty...
  31. if (lb < rshift)
  32. error = error ? error * 2 : 1;
  33. if (rshift)
  34. {
  35. BOOST_ASSERT(rshift < INT_MAX);
  36. if (bit_test(result, static_cast<unsigned>(rshift - 1)))
  37. {
  38. if (lb == rshift - 1)
  39. roundup = 1;
  40. else
  41. roundup = 2;
  42. }
  43. result >>= rshift;
  44. }
  45. if ((roundup == 2) || ((roundup == 1) && (result.backend().limbs()[0] & 1)))
  46. ++result;
  47. }
  48. return rshift;
  49. }
  50. //
  51. // Computes a^e shifted to the right so it fits in max_bits, returns how far
  52. // to the right we are shifted.
  53. //
  54. template <class I>
  55. inline I restricted_pow(cpp_int& result, const cpp_int& a, I e, I max_bits, boost::int64_t& error)
  56. {
  57. BOOST_ASSERT(&result != &a);
  58. I exp = 0;
  59. if (e == 1)
  60. {
  61. result = a;
  62. return exp;
  63. }
  64. else if (e == 2)
  65. {
  66. return restricted_multiply(result, a, a, max_bits, error);
  67. }
  68. else if (e == 3)
  69. {
  70. exp = restricted_multiply(result, a, a, max_bits, error);
  71. exp += restricted_multiply(result, result, a, max_bits, error);
  72. return exp;
  73. }
  74. I p = e / 2;
  75. exp = restricted_pow(result, a, p, max_bits, error);
  76. exp *= 2;
  77. exp += restricted_multiply(result, result, result, max_bits, error);
  78. if (e & 1)
  79. exp += restricted_multiply(result, result, a, max_bits, error);
  80. return exp;
  81. }
  82. inline int get_round_mode(const cpp_int& what, boost::int64_t location, boost::int64_t error)
  83. {
  84. //
  85. // Can we round what at /location/, if the error in what is /error/ in
  86. // units of 0.5ulp. Return:
  87. //
  88. // -1: Can't round.
  89. // 0: leave as is.
  90. // 1: tie.
  91. // 2: round up.
  92. //
  93. BOOST_ASSERT(location >= 0);
  94. BOOST_ASSERT(location < INT_MAX);
  95. boost::int64_t error_radius = error & 1 ? (1 + error) / 2 : error / 2;
  96. if (error_radius && ((int)msb(error_radius) >= location))
  97. return -1;
  98. if (bit_test(what, static_cast<unsigned>(location)))
  99. {
  100. if ((int)lsb(what) == location)
  101. return error ? -1 : 1; // Either a tie or can't round depending on whether we have any error
  102. if (!error)
  103. return 2; // no error, round up.
  104. cpp_int t = what - error_radius;
  105. if ((int)lsb(t) >= location)
  106. return -1;
  107. return 2;
  108. }
  109. else if (error)
  110. {
  111. cpp_int t = what + error_radius;
  112. return bit_test(t, static_cast<unsigned>(location)) ? -1 : 0;
  113. }
  114. return 0;
  115. }
  116. inline int get_round_mode(cpp_int& r, cpp_int& d, boost::int64_t error, const cpp_int& q)
  117. {
  118. //
  119. // Lets suppose we have an inexact division by d+delta, where the true
  120. // value for the divisor is d, and with |delta| <= error/2, then
  121. // we have calculated q and r such that:
  122. //
  123. // n r
  124. // --- = q + -----------
  125. // d + error d + error
  126. //
  127. // Rearranging for n / d we get:
  128. //
  129. // n delta*q + r
  130. // --- = q + -------------
  131. // d d
  132. //
  133. // So rounding depends on whether 2r + error * q > d.
  134. //
  135. // We return:
  136. // 0 = down down.
  137. // 1 = tie.
  138. // 2 = round up.
  139. // -1 = couldn't decide.
  140. //
  141. r <<= 1;
  142. int c = r.compare(d);
  143. if (c == 0)
  144. return error ? -1 : 1;
  145. if (c > 0)
  146. {
  147. if (error)
  148. {
  149. r -= error * q;
  150. return r.compare(d) > 0 ? 2 : -1;
  151. }
  152. return 2;
  153. }
  154. if (error)
  155. {
  156. r += error * q;
  157. return r.compare(d) < 0 ? 0 : -1;
  158. }
  159. return 0;
  160. }
  161. } // namespace cpp_bf_io_detail
  162. namespace backends {
  163. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  164. cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::operator=(const char* s)
  165. {
  166. cpp_int n;
  167. boost::intmax_t decimal_exp = 0;
  168. boost::intmax_t digits_seen = 0;
  169. static const boost::intmax_t max_digits_seen = 4 + (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count * 301L) / 1000;
  170. bool ss = false;
  171. //
  172. // Extract the sign:
  173. //
  174. if (*s == '-')
  175. {
  176. ss = true;
  177. ++s;
  178. }
  179. else if (*s == '+')
  180. ++s;
  181. //
  182. // Special cases first:
  183. //
  184. if ((std::strcmp(s, "nan") == 0) || (std::strcmp(s, "NaN") == 0) || (std::strcmp(s, "NAN") == 0))
  185. {
  186. return *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  187. }
  188. if ((std::strcmp(s, "inf") == 0) || (std::strcmp(s, "Inf") == 0) || (std::strcmp(s, "INF") == 0) || (std::strcmp(s, "infinity") == 0) || (std::strcmp(s, "Infinity") == 0) || (std::strcmp(s, "INFINITY") == 0))
  189. {
  190. *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
  191. if (ss)
  192. negate();
  193. return *this;
  194. }
  195. //
  196. // Digits before the point:
  197. //
  198. while (*s && (*s >= '0') && (*s <= '9'))
  199. {
  200. n *= 10u;
  201. n += *s - '0';
  202. if (digits_seen || (*s != '0'))
  203. ++digits_seen;
  204. ++s;
  205. }
  206. // The decimal point (we really should localise this!!)
  207. if (*s && (*s == '.'))
  208. ++s;
  209. //
  210. // Digits after the point:
  211. //
  212. while (*s && (*s >= '0') && (*s <= '9'))
  213. {
  214. n *= 10u;
  215. n += *s - '0';
  216. --decimal_exp;
  217. if (digits_seen || (*s != '0'))
  218. ++digits_seen;
  219. ++s;
  220. if (digits_seen > max_digits_seen)
  221. break;
  222. }
  223. //
  224. // Digits we're skipping:
  225. //
  226. while (*s && (*s >= '0') && (*s <= '9'))
  227. ++s;
  228. //
  229. // See if there's an exponent:
  230. //
  231. if (*s && ((*s == 'e') || (*s == 'E')))
  232. {
  233. ++s;
  234. boost::intmax_t e = 0;
  235. bool es = false;
  236. if (*s && (*s == '-'))
  237. {
  238. es = true;
  239. ++s;
  240. }
  241. else if (*s && (*s == '+'))
  242. ++s;
  243. while (*s && (*s >= '0') && (*s <= '9'))
  244. {
  245. e *= 10u;
  246. e += *s - '0';
  247. ++s;
  248. }
  249. if (es)
  250. e = -e;
  251. decimal_exp += e;
  252. }
  253. if (*s)
  254. {
  255. //
  256. // Oops unexpected input at the end of the number:
  257. //
  258. BOOST_THROW_EXCEPTION(std::runtime_error("Unable to parse string as a valid floating point number."));
  259. }
  260. if (n == 0)
  261. {
  262. // Result is necessarily zero:
  263. *this = static_cast<limb_type>(0u);
  264. return *this;
  265. }
  266. static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
  267. //
  268. // Set our working precision - this is heuristic based, we want
  269. // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations
  270. // and excessive memory usage, but we also want to avoid having to
  271. // up the computation and start again at a higher precision.
  272. // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add
  273. // one limb for good measure. This works very well for small exponents,
  274. // but for larger exponents we may may need to restart, we could add some
  275. // extra precision right from the start for larger exponents, but this
  276. // seems to be slightly slower in the *average* case:
  277. //
  278. #ifdef BOOST_MP_STRESS_IO
  279. boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32;
  280. #else
  281. boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + ((cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count % limb_bits) ? (limb_bits - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count % limb_bits) : 0) + limb_bits;
  282. #endif
  283. boost::int64_t error = 0;
  284. boost::intmax_t calc_exp = 0;
  285. boost::intmax_t final_exponent = 0;
  286. if (decimal_exp >= 0)
  287. {
  288. // Nice and simple, the result is an integer...
  289. do
  290. {
  291. cpp_int t;
  292. if (decimal_exp)
  293. {
  294. calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), decimal_exp, max_bits, error);
  295. calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(t, t, n, max_bits, error);
  296. }
  297. else
  298. t = n;
  299. final_exponent = (boost::int64_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp + calc_exp;
  300. int rshift = msb(t) - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
  301. if (rshift > 0)
  302. {
  303. final_exponent += rshift;
  304. int roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(t, rshift - 1, error);
  305. t >>= rshift;
  306. if ((roundup == 2) || ((roundup == 1) && t.backend().limbs()[0] & 1))
  307. ++t;
  308. else if (roundup < 0)
  309. {
  310. #ifdef BOOST_MP_STRESS_IO
  311. max_bits += 32;
  312. #else
  313. max_bits *= 2;
  314. #endif
  315. error = 0;
  316. continue;
  317. }
  318. }
  319. else
  320. {
  321. BOOST_ASSERT(!error);
  322. }
  323. if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
  324. {
  325. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
  326. final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
  327. }
  328. else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
  329. {
  330. // Underflow:
  331. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
  332. final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
  333. }
  334. else
  335. {
  336. exponent() = static_cast<Exponent>(final_exponent);
  337. final_exponent = 0;
  338. }
  339. copy_and_round(*this, t.backend());
  340. break;
  341. } while (true);
  342. if (ss != sign())
  343. negate();
  344. }
  345. else
  346. {
  347. // Result is the ratio of two integers: we need to organise the
  348. // division so as to produce at least an N-bit result which we can
  349. // round according to the remainder.
  350. //cpp_int d = pow(cpp_int(5), -decimal_exp);
  351. do
  352. {
  353. cpp_int d;
  354. calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -decimal_exp, max_bits, error);
  355. int shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb(n) + msb(d);
  356. final_exponent = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp - calc_exp;
  357. if (shift > 0)
  358. {
  359. n <<= shift;
  360. final_exponent -= static_cast<Exponent>(shift);
  361. }
  362. cpp_int q, r;
  363. divide_qr(n, d, q, r);
  364. int gb = msb(q);
  365. BOOST_ASSERT((gb >= static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - 1));
  366. //
  367. // Check for rounding conditions we have to
  368. // handle ourselves:
  369. //
  370. int roundup = 0;
  371. if (gb == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
  372. {
  373. // Exactly the right number of bits, use the remainder to round:
  374. roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, q);
  375. }
  376. else if (bit_test(q, gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) && ((int)lsb(q) == (gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)))
  377. {
  378. // Too many bits in q and the bits in q indicate a tie, but we can break that using r,
  379. // note that the radius of error in r is error/2 * q:
  380. int lshift = gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
  381. q >>= lshift;
  382. final_exponent += static_cast<Exponent>(lshift);
  383. BOOST_ASSERT((msb(q) >= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
  384. if (error && (r < (error / 2) * q))
  385. roundup = -1;
  386. else if (error && (r + (error / 2) * q >= d))
  387. roundup = -1;
  388. else
  389. roundup = r ? 2 : 1;
  390. }
  391. else if (error && (((error / 2) * q + r >= d) || (r < (error / 2) * q)))
  392. {
  393. // We might have been rounding up, or got the wrong quotient: can't tell!
  394. roundup = -1;
  395. }
  396. if (roundup < 0)
  397. {
  398. #ifdef BOOST_MP_STRESS_IO
  399. max_bits += 32;
  400. #else
  401. max_bits *= 2;
  402. #endif
  403. error = 0;
  404. if (shift > 0)
  405. {
  406. n >>= shift;
  407. final_exponent += static_cast<Exponent>(shift);
  408. }
  409. continue;
  410. }
  411. else if ((roundup == 2) || ((roundup == 1) && q.backend().limbs()[0] & 1))
  412. ++q;
  413. if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
  414. {
  415. // Overflow:
  416. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
  417. final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
  418. }
  419. else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
  420. {
  421. // Underflow:
  422. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
  423. final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
  424. }
  425. else
  426. {
  427. exponent() = static_cast<Exponent>(final_exponent);
  428. final_exponent = 0;
  429. }
  430. copy_and_round(*this, q.backend());
  431. if (ss != sign())
  432. negate();
  433. break;
  434. } while (true);
  435. }
  436. //
  437. // Check for scaling and/or over/under-flow:
  438. //
  439. final_exponent += exponent();
  440. if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
  441. {
  442. // Overflow:
  443. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
  444. bits() = limb_type(0);
  445. }
  446. else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
  447. {
  448. // Underflow:
  449. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  450. bits() = limb_type(0);
  451. sign() = 0;
  452. }
  453. else
  454. {
  455. exponent() = static_cast<Exponent>(final_exponent);
  456. }
  457. return *this;
  458. }
  459. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  460. std::string cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::str(std::streamsize dig, std::ios_base::fmtflags f) const
  461. {
  462. if (dig == 0)
  463. dig = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::max_digits10;
  464. bool scientific = (f & std::ios_base::scientific) == std::ios_base::scientific;
  465. bool fixed = !scientific && (f & std::ios_base::fixed);
  466. std::string s;
  467. if (exponent() <= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
  468. {
  469. // How far to left-shift in order to demormalise the mantissa:
  470. boost::intmax_t shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - (boost::intmax_t)exponent() - 1;
  471. boost::intmax_t digits_wanted = static_cast<int>(dig);
  472. boost::intmax_t base10_exp = exponent() >= 0 ? static_cast<boost::intmax_t>(std::floor(0.30103 * exponent())) : static_cast<boost::intmax_t>(std::ceil(0.30103 * exponent()));
  473. //
  474. // For fixed formatting we want /dig/ digits after the decimal point,
  475. // so if the exponent is zero, allowing for the one digit before the
  476. // decimal point, we want 1 + dig digits etc.
  477. //
  478. if (fixed)
  479. digits_wanted += 1 + base10_exp;
  480. if (scientific)
  481. digits_wanted += 1;
  482. if (digits_wanted < -1)
  483. {
  484. // Fixed precision, no significant digits, and nothing to round!
  485. s = "0";
  486. if (sign())
  487. s.insert(static_cast<std::string::size_type>(0), 1, '-');
  488. boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, true);
  489. return s;
  490. }
  491. //
  492. // power10 is the base10 exponent we need to multiply/divide by in order
  493. // to convert our denormalised number to an integer with the right number of digits:
  494. //
  495. boost::intmax_t power10 = digits_wanted - base10_exp - 1;
  496. //
  497. // If we calculate 5^power10 rather than 10^power10 we need to move
  498. // 2^power10 into /shift/
  499. //
  500. shift -= power10;
  501. cpp_int i;
  502. int roundup = 0; // 0=no rounding, 1=tie, 2=up
  503. static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
  504. //
  505. // Set our working precision - this is heuristic based, we want
  506. // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations
  507. // and excessive memory usage, but we also want to avoid having to
  508. // up the computation and start again at a higher precision.
  509. // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add
  510. // one limb for good measure. This works very well for small exponents,
  511. // but for larger exponents we add a few extra limbs to max_bits:
  512. //
  513. #ifdef BOOST_MP_STRESS_IO
  514. boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32;
  515. #else
  516. boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + ((cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count % limb_bits) ? (limb_bits - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count % limb_bits) : 0) + limb_bits;
  517. if (power10)
  518. max_bits += (msb(boost::multiprecision::detail::abs(power10)) / 8) * limb_bits;
  519. #endif
  520. do
  521. {
  522. boost::int64_t error = 0;
  523. boost::intmax_t calc_exp = 0;
  524. //
  525. // Our integer result is: bits() * 2^-shift * 5^power10
  526. //
  527. i = bits();
  528. if (shift < 0)
  529. {
  530. if (power10 >= 0)
  531. {
  532. // We go straight to the answer with all integer arithmetic,
  533. // the result is always exact and never needs rounding:
  534. BOOST_ASSERT(power10 <= (boost::intmax_t)INT_MAX);
  535. i <<= -shift;
  536. if (power10)
  537. i *= pow(cpp_int(5), static_cast<unsigned>(power10));
  538. }
  539. else if (power10 < 0)
  540. {
  541. cpp_int d;
  542. calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -power10, max_bits, error);
  543. shift += calc_exp;
  544. BOOST_ASSERT(shift < 0); // Must still be true!
  545. i <<= -shift;
  546. cpp_int r;
  547. divide_qr(i, d, i, r);
  548. roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, i);
  549. if (roundup < 0)
  550. {
  551. #ifdef BOOST_MP_STRESS_IO
  552. max_bits += 32;
  553. #else
  554. max_bits *= 2;
  555. #endif
  556. shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
  557. continue;
  558. }
  559. }
  560. }
  561. else
  562. {
  563. //
  564. // Our integer is bits() * 2^-shift * 10^power10
  565. //
  566. if (power10 > 0)
  567. {
  568. if (power10)
  569. {
  570. cpp_int t;
  571. calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), power10, max_bits, error);
  572. calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(i, i, t, max_bits, error);
  573. shift -= calc_exp;
  574. }
  575. if ((shift < 0) || ((shift == 0) && error))
  576. {
  577. // We only get here if we were asked for a crazy number of decimal digits -
  578. // more than are present in a 2^max_bits number.
  579. #ifdef BOOST_MP_STRESS_IO
  580. max_bits += 32;
  581. #else
  582. max_bits *= 2;
  583. #endif
  584. shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
  585. continue;
  586. }
  587. if (shift)
  588. {
  589. roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(i, shift - 1, error);
  590. if (roundup < 0)
  591. {
  592. #ifdef BOOST_MP_STRESS_IO
  593. max_bits += 32;
  594. #else
  595. max_bits *= 2;
  596. #endif
  597. shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
  598. continue;
  599. }
  600. i >>= shift;
  601. }
  602. }
  603. else
  604. {
  605. // We're right shifting, *and* dividing by 5^-power10,
  606. // so 5^-power10 can never be that large or we'd simply
  607. // get zero as a result, and that case is already handled above:
  608. cpp_int r;
  609. BOOST_ASSERT(-power10 < INT_MAX);
  610. cpp_int d = pow(cpp_int(5), static_cast<unsigned>(-power10));
  611. d <<= shift;
  612. divide_qr(i, d, i, r);
  613. r <<= 1;
  614. int c = r.compare(d);
  615. roundup = c < 0 ? 0 : c == 0 ? 1 : 2;
  616. }
  617. }
  618. s = i.str(0, std::ios_base::fmtflags(0));
  619. //
  620. // Check if we got the right number of digits, this
  621. // is really a test of whether we calculated the
  622. // decimal exponent correctly:
  623. //
  624. boost::intmax_t digits_got = i ? static_cast<boost::intmax_t>(s.size()) : 0;
  625. if (digits_got != digits_wanted)
  626. {
  627. base10_exp += digits_got - digits_wanted;
  628. if (fixed)
  629. digits_wanted = digits_got; // strange but true.
  630. power10 = digits_wanted - base10_exp - 1;
  631. shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
  632. if (fixed)
  633. break;
  634. roundup = 0;
  635. }
  636. else
  637. break;
  638. } while (true);
  639. //
  640. // Check whether we need to round up: note that we could equally round up
  641. // the integer /i/ above, but since we need to perform the rounding *after*
  642. // the conversion to a string and the digit count check, we might as well
  643. // do it here:
  644. //
  645. if ((roundup == 2) || ((roundup == 1) && ((s[s.size() - 1] - '0') & 1)))
  646. {
  647. boost::multiprecision::detail::round_string_up_at(s, static_cast<int>(s.size() - 1), base10_exp);
  648. }
  649. if (sign())
  650. s.insert(static_cast<std::string::size_type>(0), 1, '-');
  651. boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, false);
  652. }
  653. else
  654. {
  655. switch (exponent())
  656. {
  657. case exponent_zero:
  658. s = sign() ? "-0" : f & std::ios_base::showpos ? "+0" : "0";
  659. boost::multiprecision::detail::format_float_string(s, 0, dig, f, true);
  660. break;
  661. case exponent_nan:
  662. s = "nan";
  663. break;
  664. case exponent_infinity:
  665. s = sign() ? "-inf" : f & std::ios_base::showpos ? "+inf" : "inf";
  666. break;
  667. }
  668. }
  669. return s;
  670. }
  671. #ifdef BOOST_MSVC
  672. #pragma warning(pop)
  673. #endif
  674. } // namespace backends
  675. }} // namespace boost::multiprecision
  676. #endif