ceph/src/boost/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 http://www.boost.org/LICENSE_1_
   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{ namespace cpp_bf_io_detail{
  10
  11 #ifdef BOOST_MSVC
  12 #pragma warning(push)
  13 #pragma warning(disable:4127) // conditional expression is constant
  14 #endif
  15
  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
 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       }
 353       while(true);
 354
 355       if(ss != sign())
 356          negate();
 357    }
 358    else
 359    {
 360       // Result is the ratio of two integers: we need to organise the
 361       // division so as to produce at least an N-bit result which we can
 362       // round according to the remainder.
 363       //cpp_int d = pow(cpp_int(5), -decimal_exp);
 364       do
 365       {
 366          cpp_int d;
 367          calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -decimal_exp, max_bits, error);
 368          int shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb(n) + msb(d);
 369          final_exponent = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp - calc_exp;
 370          if(shift > 0)
 371          {
 372             n <<= shift;
 373             final_exponent -= static_cast<Exponent>(shift);
 374          }
 375          cpp_int q, r;
 376          divide_qr(n, d, q, r);
 377          int gb = msb(q);
 378          BOOST_ASSERT((gb >= static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - 1));
 379          //
 380          // Check for rounding conditions we have to
 381          // handle ourselves:
 382          //
 383          int roundup = 0;
 384          if(gb == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
 385          {
 386             // Exactly the right number of bits, use the remainder to round:
 387             roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, q);
 388          }
 389          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)))
 390          {
 391             // Too many bits in q and the bits in q indicate a tie, but we can break that using r,
 392             // note that the radius of error in r is error/2 * q:
 393             int lshift = gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
 394             q >>= lshift;
 395             final_exponent += static_cast<Exponent>(lshift);
 396             BOOST_ASSERT((msb(q) >= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
 397             if(error && (r < (error / 2) * q))
 398                roundup = -1;
 399             else if(error && (r + (error / 2) * q >= d))
 400                roundup = -1;
 401             else
 402                roundup = r ? 2 : 1;
 403          }
 404          else if(error && (((error / 2) * q + r >= d) || (r < (error / 2) * q)))
 405          {
 406             // We might have been rounding up, or got the wrong quotient: can't tell!
 407             roundup = -1;
 408          }
 409          if(roundup < 0)
 410          {
 411 #ifdef BOOST_MP_STRESS_IO
 412             max_bits += 32;
 413 #else
 414             max_bits *= 2;
 415 #endif
 416             error = 0;
 417             if(shift > 0)
 418             {
 419                n >>= shift;
 420                final_exponent += static_cast<Exponent>(shift);
 421             }
 422             continue;
 423          }
 424          else if((roundup == 2) || ((roundup == 1) && q.backend().limbs()[0] & 1))
 425             ++q;
 426          if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
 427          {
 428             // Overflow:
 429             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
 430             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
 431          }
 432          else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
 433          {
 434             // Underflow:
 435             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
 436             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
 437          }
 438          else
 439          {
 440             exponent() = static_cast<Exponent>(final_exponent);
 441             final_exponent = 0;
 442          }
 443          copy_and_round(*this, q.backend());
 444          if(ss != sign())
 445             negate();
 446          break;
 447       }
 448       while(true);
 449    }
 450    //
 451    // Check for scaling and/or over/under-flow:
 452    //
 453    final_exponent += exponent();
 454    if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
 455    {
 456       // Overflow:
 457       exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
 458       bits() = limb_type(0);
 459    }
 460    else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
 461    {
 462       // Underflow:
 463       exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
 464       bits() = limb_type(0);
 465       sign() = 0;
 466    }
 467    else
 468    {
 469       exponent() = static_cast<Exponent>(final_exponent);
 470    }
 471    return *this;
 472 }
 473
 474 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 475 std::string cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::str(std::streamsize dig, std::ios_base::fmtflags f) const
 476 {
 477    if(dig == 0)
 478       dig = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::max_digits10;
 479
 480    bool scientific = (f & std::ios_base::scientific) == std::ios_base::scientific;
 481    bool fixed = !scientific && (f & std::ios_base::fixed);
 482
 483    std::string s;
 484
 485    if(exponent() <= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
 486    {
 487       // How far to left-shift in order to demormalise the mantissa:
 488       boost::intmax_t shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1;
 489       boost::intmax_t digits_wanted = static_cast<int>(dig);
 490       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()));
 491       //
 492       // For fixed formatting we want /dig/ digits after the decimal point,
 493       // so if the exponent is zero, allowing for the one digit before the
 494       // decimal point, we want 1 + dig digits etc.
 495       //
 496       if(fixed)
 497          digits_wanted += 1 + base10_exp;
 498       if(scientific)
 499          digits_wanted += 1;
 500       if(digits_wanted < -1)
 501       {
 502          // Fixed precision, no significant digits, and nothing to round!
 503          s = "0";
 504          if(sign())
 505             s.insert(static_cast<std::string::size_type>(0), 1, '-');
 506          boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, true);
 507          return s;
 508       }
 509       //
 510       // power10 is the base10 exponent we need to multiply/divide by in order
 511       // to convert our denormalised number to an integer with the right number of digits:
 512       //
 513       boost::intmax_t power10 = digits_wanted - base10_exp - 1;
 514       //
 515       // If we calculate 5^power10 rather than 10^power10 we need to move
 516       // 2^power10 into /shift/
 517       //
 518       shift -= power10;
 519       cpp_int i;
 520       int roundup = 0; // 0=no rounding, 1=tie, 2=up
 521       static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
 522       //
 523       // Set our working precision - this is heuristic based, we want
 524       // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations
 525       // and excessive memory usage, but we also want to avoid having to
 526       // up the computation and start again at a higher precision.
 527       // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add
 528       // one limb for good measure.  This works very well for small exponents,
 529       // but for larger exponents we add a few extra limbs to max_bits:
 530       //
 531 #ifdef BOOST_MP_STRESS_IO
 532       boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32;
 533 #else
 534       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;
 535       if(power10)
 536          max_bits += (msb(boost::multiprecision::detail::abs(power10)) / 8) * limb_bits;
 537 #endif
 538       do
 539       {
 540          boost::int64_t error = 0;
 541          boost::intmax_t calc_exp = 0;
 542          //
 543          // Our integer result is: bits() * 2^-shift * 5^power10
 544          //
 545          i = bits();
 546          if(shift < 0)
 547          {
 548             if(power10 >= 0)
 549             {
 550                // We go straight to the answer with all integer arithmetic,
 551                // the result is always exact and never needs rounding:
 552                BOOST_ASSERT(power10 <= (boost::intmax_t)INT_MAX);
 553                i <<= -shift;
 554                if(power10)
 555                   i *= pow(cpp_int(5), static_cast<unsigned>(power10));
 556             }
 557             else if(power10 < 0)
 558             {
 559                cpp_int d;
 560                calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -power10, max_bits, error);
 561                shift += calc_exp;
 562                BOOST_ASSERT(shift < 0); // Must still be true!
 563                i <<= -shift;
 564                cpp_int r;
 565                divide_qr(i, d, i, r);
 566                roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, i);
 567                if(roundup < 0)
 568                {
 569 #ifdef BOOST_MP_STRESS_IO
 570                   max_bits += 32;
 571 #else
 572                   max_bits *= 2;
 573 #endif
 574                   shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
 575                   continue;
 576                }
 577             }
 578          }
 579          else
 580          {
 581             //
 582             // Our integer is bits() * 2^-shift * 10^power10
 583             //
 584             if(power10 > 0)
 585             {
 586                if(power10)
 587                {
 588                   cpp_int t;
 589                   calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), power10, max_bits, error);
 590                   calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(i, i, t, max_bits, error);
 591                   shift -= calc_exp;
 592                }
 593                if((shift < 0) || ((shift == 0) && error))
 594                {
 595                   // We only get here if we were asked for a crazy number of decimal digits -
 596                   // more than are present in a 2^max_bits number.
 597 #ifdef BOOST_MP_STRESS_IO
 598                   max_bits += 32;
 599 #else
 600                   max_bits *= 2;
 601 #endif
 602                   shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
 603                   continue;
 604                }
 605                if(shift)
 606                {
 607                   roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(i, shift - 1, error);
 608                   if(roundup < 0)
 609                   {
 610 #ifdef BOOST_MP_STRESS_IO
 611                      max_bits += 32;
 612 #else
 613                      max_bits *= 2;
 614 #endif
 615                      shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
 616                      continue;
 617                   }
 618                   i >>= shift;
 619                }
 620             }
 621             else
 622             {
 623                // We're right shifting, *and* dividing by 5^-power10,
 624                // so 5^-power10 can never be that large or we'd simply
 625                // get zero as a result, and that case is already handled above:
 626                cpp_int r;
 627                BOOST_ASSERT(-power10 < INT_MAX);
 628                cpp_int d = pow(cpp_int(5), static_cast<unsigned>(-power10));
 629                d <<= shift;
 630                divide_qr(i, d, i, r);
 631                r <<= 1;
 632                int c = r.compare(d);
 633                roundup = c < 0 ? 0 : c == 0 ? 1 : 2;
 634             }
 635          }
 636          s = i.str(0, std::ios_base::fmtflags(0));
 637          //
 638          // Check if we got the right number of digits, this
 639          // is really a test of whether we calculated the
 640          // decimal exponent correctly:
 641          //
 642          boost::intmax_t digits_got = i ? static_cast<boost::intmax_t>(s.size()) : 0;
 643          if(digits_got != digits_wanted)
 644          {
 645             base10_exp += digits_got - digits_wanted;
 646             if(fixed)
 647                digits_wanted = digits_got;  // strange but true.
 648             power10 = digits_wanted - base10_exp - 1;
 649             shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
 650             if(fixed)
 651                break;
 652             roundup = 0;
 653          }
 654          else
 655             break;
 656       }
 657       while(true);
 658       //
 659       // Check whether we need to round up: note that we could equally round up
 660       // the integer /i/ above, but since we need to perform the rounding *after*
 661       // the conversion to a string and the digit count check, we might as well
 662       // do it here:
 663       //
 664       if((roundup == 2) || ((roundup == 1) && ((s[s.size() - 1] - '0') & 1)))
 665       {
 666          boost::multiprecision::detail::round_string_up_at(s, static_cast<int>(s.size() - 1), base10_exp);
 667       }
 668
 669       if(sign())
 670          s.insert(static_cast<std::string::size_type>(0), 1, '-');
 671
 672       boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, false);
 673    }
 674    else
 675    {
 676       switch(exponent())
 677       {
 678       case exponent_zero:
 679          s = sign() ? "-0" : f & std::ios_base::showpos ? "+0" : "0";
 680          boost::multiprecision::detail::format_float_string(s, 0, dig, f, true);
 681          break;
 682       case exponent_nan:
 683          s = "nan";
 684          break;
 685       case exponent_infinity:
 686          s = sign() ? "-inf" : f & std::ios_base::showpos ? "+inf" : "inf";
 687          break;
 688       }
 689    }
 690    return s;
 691 }
 692
 693 #ifdef BOOST_MSVC
 694 #pragma warning(pop)
 695 #endif
 696
 697 }}} // namespaces
 698
 699 #endif
 700