boost/multiprecision/cpp_bin_float/io.hpp

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