io.hpp 25 KB

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