multiprecision/cpp_bin_float/io.hpp

0001 ///////////////////////////////////////////////////////////////
0002 //  Copyright 2013 John Maddock. Distributed under the Boost
0003 //  Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
0005
0006 #ifndef BOOST_MP_CPP_BIN_FLOAT_IO_HPP
0007 #define BOOST_MP_CPP_BIN_FLOAT_IO_HPP
0008
0009 #include <boost/multiprecision/detail/no_exceptions_support.hpp>
0010 #include <boost/multiprecision/detail/assert.hpp>
0011
0012 namespace boost { namespace multiprecision {
0013 namespace cpp_bf_io_detail {
0014
0015 //
0016 // Multiplies a by b and shifts the result so it fits inside max_bits bits,
0017 // returns by how much the result was shifted.
0018 //
0019 template <class I>
0020 inline I restricted_multiply(cpp_int& result, const cpp_int& a, const cpp_int& b, I max_bits, std::int64_t& error)
0021 {
0022    using local_integral_type = I;
0023
0024    result   = a * b;
0025    local_integral_type gb     = static_cast<local_integral_type>(msb(result));
0026    local_integral_type rshift = 0;
0027    if (gb > max_bits)
0028    {
0029       rshift = gb - max_bits;
0030       local_integral_type lb      = static_cast<local_integral_type>(lsb(result));
0031       int roundup = 0;
0032       // The error rate increases by the error of both a and b,
0033       // this may be overly pessimistic in many case as we're assuming
0034       // that a and b have the same level of uncertainty...
0035       if (lb < rshift)
0036          error = error ? error * 2 : 1;
0037       if (rshift)
0038       {
0039          BOOST_MP_ASSERT(rshift < INT_MAX);
0040          if (bit_test(result, static_cast<unsigned>(rshift - 1)))
0041          {
0042             if (lb == rshift - 1)
0043                roundup = 1;
0044             else
0045                roundup = 2;
0046          }
0047          result >>= rshift;
0048       }
0049       if ((roundup == 2) || ((roundup == 1) && (result.backend().limbs()[0] & 1)))
0050          ++result;
0051    }
0052    return rshift;
0053 }
0054 //
0055 // Computes a^e shifted to the right so it fits in max_bits, returns how far
0056 // to the right we are shifted.
0057 //
0058 template <class I>
0059 inline I restricted_pow(cpp_int& result, const cpp_int& a, I e, I max_bits, std::int64_t& error)
0060 {
0061    BOOST_MP_ASSERT(&result != &a);
0062    I exp = 0;
0063    if (e == 1)
0064    {
0065       result = a;
0066       return exp;
0067    }
0068    else if (e == 2)
0069    {
0070       return restricted_multiply(result, a, a, max_bits, error);
0071    }
0072    else if (e == 3)
0073    {
0074       exp = restricted_multiply(result, a, a, max_bits, error);
0075       exp += restricted_multiply(result, result, a, max_bits, error);
0076       return exp;
0077    }
0078    I p = e / 2;
0079    exp = restricted_pow(result, a, p, max_bits, error);
0080    exp *= 2;
0081    exp += restricted_multiply(result, result, result, max_bits, error);
0082    if (e & 1)
0083       exp += restricted_multiply(result, result, a, max_bits, error);
0084    return exp;
0085 }
0086
0087 inline int get_round_mode(const cpp_int& what, std::int64_t location, std::int64_t error)
0088 {
0089    //
0090    // Can we round what at /location/, if the error in what is /error/ in
0091    // units of 0.5ulp.  Return:
0092    //
0093    // -1: Can't round.
0094    //  0: leave as is.
0095    //  1: tie.
0096    //  2: round up.
0097    //
0098    BOOST_MP_ASSERT(location >= 0);
0099    BOOST_MP_ASSERT(location < INT_MAX);
0100    std::int64_t error_radius = error & 1 ? (1 + error) / 2 : error / 2;
0101    if (error_radius && (static_cast<int>(msb(error_radius)) >= location))
0102       return -1;
0103    if (bit_test(what, static_cast<unsigned>(location)))
0104    {
0105       if (static_cast<int>(lsb(what)) == location)
0106          return error ? -1 : 1; // Either a tie or can't round depending on whether we have any error
0107       if (!error)
0108          return 2; // no error, round up.
0109       cpp_int t = what - error_radius;
0110       if (static_cast<int>(lsb(t)) >= location)
0111          return -1;
0112       return 2;
0113    }
0114    else if (error)
0115    {
0116       cpp_int t = what + error_radius;
0117       return bit_test(t, static_cast<unsigned>(location)) ? -1 : 0;
0118    }
0119    return 0;
0120 }
0121
0122 inline int get_round_mode(cpp_int& r, cpp_int& d, std::int64_t error, const cpp_int& q)
0123 {
0124    //
0125    // Lets suppose we have an inexact division by d+delta, where the true
0126    // value for the divisor is d, and with |delta| <= error/2, then
0127    // we have calculated q and r such that:
0128    //
0129    // n                  r
0130    // ---       = q + -----------
0131    // d + error        d + error
0132    //
0133    // Rearranging for n / d we get:
0134    //
0135    //    n         delta*q + r
0136    //   --- = q + -------------
0137    //    d              d
0138    //
0139    // So rounding depends on whether 2r + error * q > d.
0140    //
0141    // We return:
0142    //  0 = down down.
0143    //  1 = tie.
0144    //  2 = round up.
0145    // -1 = couldn't decide.
0146    //
0147    r <<= 1;
0148    int c = r.compare(d);
0149    if (c == 0)
0150       return error ? -1 : 1;
0151    if (c > 0)
0152    {
0153       if (error)
0154       {
0155          r -= error * q;
0156          return r.compare(d) > 0 ? 2 : -1;
0157       }
0158       return 2;
0159    }
0160    if (error)
0161    {
0162       r += error * q;
0163       return r.compare(d) < 0 ? 0 : -1;
0164    }
0165    return 0;
0166 }
0167
0168 } // namespace cpp_bf_io_detail
0169
0170 namespace backends {
0171
0172 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
0173 cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::operator=(const char* s)
0174 {
0175              cpp_int       n;
0176              std::intmax_t decimal_exp     = 0;
0177              std::intmax_t digits_seen     = 0;
0178    constexpr std::intmax_t max_digits_seen = 4 + (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count * 301L) / 1000;
0179              bool          ss              = false;
0180    //
0181    // Extract the sign:
0182    //
0183    if (*s == '-')
0184    {
0185       ss = true;
0186       ++s;
0187    }
0188    else if (*s == '+')
0189       ++s;
0190    //
0191    // Special cases first:
0192    //
0193    if ((std::strcmp(s, "nan") == 0) || (std::strcmp(s, "NaN") == 0) || (std::strcmp(s, "NAN") == 0))
0194    {
0195       return *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
0196    }
0197    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))
0198    {
0199       *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
0200       if (ss)
0201          negate();
0202       return *this;
0203    }
0204    //
0205    // Digits before the point:
0206    //
0207    while (*s && (*s >= '0') && (*s <= '9'))
0208    {
0209       n *= 10u;
0210       n += *s - '0';
0211       if (digits_seen || (*s != '0'))
0212          ++digits_seen;
0213       ++s;
0214    }
0215    // The decimal point (we really should localise this!!)
0216    if (*s && (*s == '.'))
0217       ++s;
0218    //
0219    // Digits after the point:
0220    //
0221    while (*s && (*s >= '0') && (*s <= '9'))
0222    {
0223       n *= 10u;
0224       n += *s - '0';
0225       --decimal_exp;
0226       if (digits_seen || (*s != '0'))
0227          ++digits_seen;
0228       ++s;
0229       if (digits_seen > max_digits_seen)
0230          break;
0231    }
0232    //
0233    // Digits we're skipping:
0234    //
0235    while (*s && (*s >= '0') && (*s <= '9'))
0236       ++s;
0237    //
0238    // See if there's an exponent:
0239    //
0240    if (*s && ((*s == 'e') || (*s == 'E')))
0241    {
0242       ++s;
0243       std::intmax_t e  = 0;
0244       bool            es = false;
0245       if (*s && (*s == '-'))
0246       {
0247          es = true;
0248          ++s;
0249       }
0250       else if (*s && (*s == '+'))
0251          ++s;
0252       while (*s && (*s >= '0') && (*s <= '9'))
0253       {
0254          e *= 10u;
0255          e += *s - '0';
0256          ++s;
0257       }
0258       if (es)
0259          e = -e;
0260       decimal_exp += e;
0261    }
0262    if (*s)
0263    {
0264       //
0265       // Oops unexpected input at the end of the number:
0266       //
0267       BOOST_MP_THROW_EXCEPTION(std::runtime_error("Unable to parse string as a valid floating point number."));
0268    }
0269    if (n == 0)
0270    {
0271       // Result is necessarily zero:
0272       *this = static_cast<limb_type>(0u);
0273       return *this;
0274    }
0275
0276    constexpr std::size_t limb_bits = sizeof(limb_type) * CHAR_BIT;
0277    //
0278    // Set our working precision - this is heuristic based, we want
0279    // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations
0280    // and excessive memory usage, but we also want to avoid having to
0281    // up the computation and start again at a higher precision.
0282    // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add
0283    // one limb for good measure.  This works very well for small exponents,
0284    // but for larger exponents we may may need to restart, we could add some
0285    // extra precision right from the start for larger exponents, but this
0286    // seems to be slightly slower in the *average* case:
0287    //
0288 #ifdef BOOST_MP_STRESS_IO
0289    std::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32;
0290 #else
0291    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;
0292 #endif
0293    std::int64_t  error          = 0;
0294    std::intmax_t calc_exp       = 0;
0295    std::intmax_t final_exponent = 0;
0296
0297    if (decimal_exp >= 0)
0298    {
0299       // Nice and simple, the result is an integer...
0300       do
0301       {
0302          cpp_int t;
0303          if (decimal_exp)
0304          {
0305             calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), decimal_exp, max_bits, error);
0306             calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(t, t, n, max_bits, error);
0307          }
0308          else
0309             t = n;
0310          final_exponent = (std::int64_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp + calc_exp;
0311          std::ptrdiff_t rshift     = static_cast<std::ptrdiff_t>(static_cast<std::ptrdiff_t>(msb(t)) - static_cast<std::ptrdiff_t>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) + 1);
0312          if (rshift > 0)
0313          {
0314             final_exponent += rshift;
0315             int roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(t, rshift - 1, error);
0316             t >>= rshift;
0317             if ((roundup == 2) || ((roundup == 1) && t.backend().limbs()[0] & 1))
0318                ++t;
0319             else if (roundup < 0)
0320             {
0321 #ifdef BOOST_MP_STRESS_IO
0322                max_bits += 32;
0323 #else
0324                max_bits *= 2;
0325 #endif
0326                error = 0;
0327                continue;
0328             }
0329          }
0330          else
0331          {
0332             BOOST_MP_ASSERT(!error);
0333          }
0334          if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
0335          {
0336             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
0337             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
0338          }
0339          else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
0340          {
0341             // Underflow:
0342             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
0343             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
0344          }
0345          else
0346          {
0347             exponent()     = static_cast<Exponent>(final_exponent);
0348             final_exponent = 0;
0349          }
0350          copy_and_round(*this, t.backend());
0351          break;
0352       } while (true);
0353
0354       if (ss != sign())
0355          negate();
0356    }
0357    else
0358    {
0359       // Result is the ratio of two integers: we need to organise the
0360       // division so as to produce at least an N-bit result which we can
0361       // round according to the remainder.
0362       //cpp_int d = pow(cpp_int(5), -decimal_exp);
0363       do
0364       {
0365          cpp_int d;
0366          calc_exp                   = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -decimal_exp, max_bits, error);
0367          const std::ptrdiff_t shift = static_cast<std::ptrdiff_t>(static_cast<std::ptrdiff_t>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - static_cast<std::ptrdiff_t>(msb(n)) + static_cast<std::ptrdiff_t>(msb(d)));
0368          final_exponent             = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp - calc_exp;
0369          if (shift > 0)
0370          {
0371             n <<= shift;
0372             final_exponent -= static_cast<Exponent>(shift);
0373          }
0374          cpp_int q, r;
0375          divide_qr(n, d, q, r);
0376          std::ptrdiff_t gb = static_cast<std::ptrdiff_t>(msb(q));
0377          BOOST_MP_ASSERT((gb >= static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - 1));
0378          //
0379          // Check for rounding conditions we have to
0380          // handle ourselves:
0381          //
0382          int roundup = 0;
0383          if (gb == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
0384          {
0385             // Exactly the right number of bits, use the remainder to round:
0386             roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, q);
0387          }
0388          else if (bit_test(q, static_cast<std::size_t>(gb - static_cast<std::ptrdiff_t>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))) && (static_cast<std::ptrdiff_t>(lsb(q)) == static_cast<std::ptrdiff_t>(gb - static_cast<std::ptrdiff_t>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))))
0389          {
0390             // Too many bits in q and the bits in q indicate a tie, but we can break that using r,
0391             // note that the radius of error in r is error/2 * q:
0392             std::ptrdiff_t lshift = gb - static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) + 1;
0393             q >>= lshift;
0394             final_exponent += static_cast<Exponent>(lshift);
0395             BOOST_MP_ASSERT((msb(q) >= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
0396             if (error && (r < (error / 2) * q))
0397                roundup = -1;
0398             else if (error && (r + (error / 2) * q >= d))
0399                roundup = -1;
0400             else
0401                roundup = r ? 2 : 1;
0402          }
0403          else if (error && (((error / 2) * q + r >= d) || (r < (error / 2) * q)))
0404          {
0405             // We might have been rounding up, or got the wrong quotient: can't tell!
0406             roundup = -1;
0407          }
0408          if (roundup < 0)
0409          {
0410 #ifdef BOOST_MP_STRESS_IO
0411             max_bits += 32;
0412 #else
0413             max_bits *= 2;
0414 #endif
0415             error = 0;
0416             if (shift > 0)
0417             {
0418                n >>= shift;
0419                final_exponent += static_cast<Exponent>(shift);
0420             }
0421             continue;
0422          }
0423          else if ((roundup == 2) || ((roundup == 1) && q.backend().limbs()[0] & 1))
0424             ++q;
0425          if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
0426          {
0427             // Overflow:
0428             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
0429             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
0430          }
0431          else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
0432          {
0433             // Underflow:
0434             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
0435             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
0436          }
0437          else
0438          {
0439             exponent()     = static_cast<Exponent>(final_exponent);
0440             final_exponent = 0;
0441          }
0442          copy_and_round(*this, q.backend());
0443          if (ss != sign())
0444             negate();
0445          break;
0446       } while (true);
0447    }
0448    //
0449    // Check for scaling and/or over/under-flow:
0450    //
0451    final_exponent += exponent();
0452    if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
0453    {
0454       // Overflow:
0455       exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
0456       bits()     = limb_type(0);
0457    }
0458    else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
0459    {
0460       // Underflow:
0461       exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
0462       bits()     = limb_type(0);
0463       sign()     = 0;
0464    }
0465    else
0466    {
0467       exponent() = static_cast<Exponent>(final_exponent);
0468    }
0469    return *this;
0470 }
0471
0472 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
0473 std::string cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::str(std::streamsize dig, std::ios_base::fmtflags f) const
0474 {
0475    bool scientific = (f & std::ios_base::scientific) == std::ios_base::scientific;
0476    bool fixed      = !scientific && (f & std::ios_base::fixed);
0477
0478    if (dig == 0 && !fixed)
0479       dig = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::max_digits10;
0480
0481    std::string s;
0482
0483    if (exponent() <= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
0484    {
0485       // How far to left-shift in order to demormalise the mantissa:
0486       std::intmax_t shift         = (std::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - (std::intmax_t)exponent() - 1;
0487       std::intmax_t digits_wanted = static_cast<int>(dig);
0488       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()));
0489       //
0490       // For fixed formatting we want /dig/ digits after the decimal point,
0491       // so if the exponent is zero, allowing for the one digit before the
0492       // decimal point, we want 1 + dig digits etc.
0493       //
0494       if (fixed)
0495          digits_wanted += 1 + base10_exp;
0496       if (scientific)
0497          digits_wanted += 1;
0498       if (digits_wanted < -1)
0499       {
0500          // Fixed precision, no significant digits, and nothing to round!
0501          s = "0";
0502          if (sign())
0503             s.insert(static_cast<std::string::size_type>(0), 1, '-');
0504          boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, true);
0505          return s;
0506       }
0507       //
0508       // power10 is the base10 exponent we need to multiply/divide by in order
0509       // to convert our denormalised number to an integer with the right number of digits:
0510       //
0511       std::intmax_t power10 = digits_wanted - base10_exp - 1;
0512       //
0513       // If we calculate 5^power10 rather than 10^power10 we need to move
0514       // 2^power10 into /shift/
0515       //
0516       shift -= power10;
0517       cpp_int               i;
0518       int                   roundup   = 0; // 0=no rounding, 1=tie, 2=up
0519       constexpr std::size_t limb_bits = sizeof(limb_type) * CHAR_BIT;
0520       //
0521       // Set our working precision - this is heuristic based, we want
0522       // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations
0523       // and excessive memory usage, but we also want to avoid having to
0524       // up the computation and start again at a higher precision.
0525       // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add
0526       // one limb for good measure.  This works very well for small exponents,
0527       // but for larger exponents we add a few extra limbs to max_bits:
0528       //
0529 #ifdef BOOST_MP_STRESS_IO
0530       std::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32;
0531 #else
0532       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;
0533       if (power10)
0534       {
0535          const std::uintmax_t abs_power10 = static_cast<std::uintmax_t>(boost::multiprecision::detail::abs(power10));
0536          const std::intmax_t  msb_div8    = static_cast<std::intmax_t>(msb(abs_power10) / 8u);
0537
0538          max_bits += (msb_div8 * static_cast<std::intmax_t>(limb_bits));
0539       }
0540 #endif
0541       do
0542       {
0543          std::int64_t  error    = 0;
0544          std::intmax_t calc_exp = 0;
0545          //
0546          // Our integer result is: bits() * 2^-shift * 5^power10
0547          //
0548          i = bits();
0549          if (shift < 0)
0550          {
0551             if (power10 >= 0)
0552             {
0553                // We go straight to the answer with all integer arithmetic,
0554                // the result is always exact and never needs rounding:
0555                BOOST_MP_ASSERT(power10 <= (std::intmax_t)INT_MAX);
0556                i <<= -shift;
0557                if (power10)
0558                   i *= pow(cpp_int(5), static_cast<unsigned>(power10));
0559             }
0560             else if (power10 < 0)
0561             {
0562                cpp_int d;
0563                calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -power10, max_bits, error);
0564                shift += calc_exp;
0565                BOOST_MP_ASSERT(shift < 0); // Must still be true!
0566                i <<= -shift;
0567                cpp_int r;
0568                divide_qr(i, d, i, r);
0569                roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, i);
0570                if (roundup < 0)
0571                {
0572 #ifdef BOOST_MP_STRESS_IO
0573                   max_bits += 32;
0574 #else
0575                   max_bits *= 2;
0576 #endif
0577                   shift = (std::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
0578                   continue;
0579                }
0580             }
0581          }
0582          else
0583          {
0584             //
0585             // Our integer is bits() * 2^-shift * 10^power10
0586             //
0587             if (power10 > 0)
0588             {
0589                if (power10)
0590                {
0591                   cpp_int t;
0592                   calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), power10, max_bits, error);
0593                   calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(i, i, t, max_bits, error);
0594                   shift -= calc_exp;
0595                }
0596                if ((shift < 0) || ((shift == 0) && error))
0597                {
0598                   // We only get here if we were asked for a crazy number of decimal digits -
0599                   // more than are present in a 2^max_bits number.
0600 #ifdef BOOST_MP_STRESS_IO
0601                   max_bits += 32;
0602 #else
0603                   max_bits *= 2;
0604 #endif
0605                   shift = (std::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
0606                   continue;
0607                }
0608                if (shift)
0609                {
0610                   roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(i, shift - 1, error);
0611                   if (roundup < 0)
0612                   {
0613 #ifdef BOOST_MP_STRESS_IO
0614                      max_bits += 32;
0615 #else
0616                      max_bits *= 2;
0617 #endif
0618                      shift = (std::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
0619                      continue;
0620                   }
0621                   i >>= shift;
0622                }
0623             }
0624             else
0625             {
0626                // We're right shifting, *and* dividing by 5^-power10,
0627                // so 5^-power10 can never be that large or we'd simply
0628                // get zero as a result, and that case is already handled above:
0629                cpp_int r;
0630                BOOST_MP_ASSERT(-power10 < INT_MAX);
0631                cpp_int d = pow(cpp_int(5), static_cast<unsigned>(-power10));
0632                d <<= shift;
0633                divide_qr(i, d, i, r);
0634                r <<= 1;
0635                int c   = r.compare(d);
0636                roundup = c < 0 ? 0 : c == 0 ? 1 : 2;
0637             }
0638          }
0639          s = i.str(0, std::ios_base::fmtflags(0));
0640          //
0641          // Check if we got the right number of digits, this
0642          // is really a test of whether we calculated the
0643          // decimal exponent correctly:
0644          //
0645          std::intmax_t digits_got = i ? static_cast<std::intmax_t>(s.size()) : 0;
0646          if (digits_got != digits_wanted)
0647          {
0648             base10_exp += digits_got - digits_wanted;
0649             if (fixed)
0650                digits_wanted = digits_got; // strange but true.
0651             power10 = digits_wanted - base10_exp - 1;
0652             shift   = (std::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
0653             if (fixed)
0654                break;
0655             roundup = 0;
0656          }
0657          else
0658             break;
0659       } while (true);
0660       //
0661       // Check whether we need to round up: note that we could equally round up
0662       // the integer /i/ above, but since we need to perform the rounding *after*
0663       // the conversion to a string and the digit count check, we might as well
0664       // do it here:
0665       //
0666       if ((roundup == 2) || ((roundup == 1) && ((s[s.size() - 1] - '0') & 1)))
0667       {
0668          boost::multiprecision::detail::round_string_up_at(s, static_cast<int>(s.size() - 1), base10_exp);
0669       }
0670
0671       if (sign())
0672          s.insert(static_cast<std::string::size_type>(0), 1, '-');
0673
0674       boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, false);
0675    }
0676    else
0677    {
0678       switch (exponent())
0679       {
0680       case exponent_zero:
0681          s = sign() ? "-0" : f & std::ios_base::showpos ? "+0" : "0";
0682          boost::multiprecision::detail::format_float_string(s, 0, dig, f, true);
0683          break;
0684       case exponent_nan:
0685          s = "nan";
0686          break;
0687       case exponent_infinity:
0688          s = sign() ? "-inf" : f & std::ios_base::showpos ? "+inf" : "inf";
0689          break;
0690       }
0691    }
0692    return s;
0693 }
0694
0695 } // namespace backends
0696 }} // namespace boost::multiprecision
0697
0698 #endif