Chris@102: ///////////////////////////////////////////////////////////////
Chris@102: //  Copyright 2013 John Maddock. Distributed under the Boost
Chris@102: //  Software License, Version 1.0. (See accompanying file
Chris@102: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
Chris@102: 
Chris@102: #ifndef BOOST_MP_CPP_BIN_FLOAT_IO_HPP
Chris@102: #define BOOST_MP_CPP_BIN_FLOAT_IO_HPP
Chris@102: 
Chris@102: namespace boost{ namespace multiprecision{ namespace cpp_bf_io_detail{
Chris@102: 
Chris@102: //
Chris@102: // Multiplies a by b and shifts the result so it fits inside max_bits bits, 
Chris@102: // returns by how much the result was shifted.
Chris@102: //
Chris@102: template <class I>
Chris@102: inline I restricted_multiply(cpp_int& result, const cpp_int& a, const cpp_int& b, I max_bits, boost::int64_t& error)
Chris@102: {
Chris@102:    result = a * b;
Chris@102:    I gb = msb(result);
Chris@102:    I rshift = 0;
Chris@102:    if(gb > max_bits)
Chris@102:    {
Chris@102:       rshift = gb - max_bits;
Chris@102:       I lb = lsb(result);
Chris@102:       int roundup = 0;
Chris@102:       // The error rate increases by the error of both a and b, 
Chris@102:       // this may be overly pessimistic in many case as we're assuming
Chris@102:       // that a and b have the same level of uncertainty...
Chris@102:       if(lb < rshift)
Chris@102:          error = error ? error * 2 : 1;
Chris@102:       if(rshift)
Chris@102:       {
Chris@102:          BOOST_ASSERT(rshift < INT_MAX);
Chris@102:          if(bit_test(result, static_cast<unsigned>(rshift - 1)))
Chris@102:          {
Chris@102:             if(lb == rshift - 1)
Chris@102:                roundup = 1;
Chris@102:             else
Chris@102:                roundup = 2;
Chris@102:          }
Chris@102:          result >>= rshift;
Chris@102:       }
Chris@102:       if((roundup == 2) || ((roundup == 1) && (result.backend().limbs()[0] & 1)))
Chris@102:          ++result;
Chris@102:    }
Chris@102:    return rshift;
Chris@102: }
Chris@102: //
Chris@102: // Computes a^e shifted to the right so it fits in max_bits, returns how far
Chris@102: // to the right we are shifted.
Chris@102: //
Chris@102: template <class I>
Chris@102: inline I restricted_pow(cpp_int& result, const cpp_int& a, I e, I max_bits, boost::int64_t& error)
Chris@102: {
Chris@102:    BOOST_ASSERT(&result != &a);
Chris@102:    I exp = 0;
Chris@102:    if(e == 1)
Chris@102:    {
Chris@102:       result = a;
Chris@102:       return exp;
Chris@102:    }
Chris@102:    else if(e == 2)
Chris@102:    {
Chris@102:       return restricted_multiply(result, a, a, max_bits, error);
Chris@102:    }
Chris@102:    else if(e == 3)
Chris@102:    {
Chris@102:       exp = restricted_multiply(result, a, a, max_bits, error);
Chris@102:       exp += restricted_multiply(result, result, a, max_bits, error);
Chris@102:       return exp;
Chris@102:    }
Chris@102:    I p = e / 2;
Chris@102:    exp = restricted_pow(result, a, p, max_bits, error);
Chris@102:    exp *= 2;
Chris@102:    exp += restricted_multiply(result, result, result, max_bits, error);
Chris@102:    if(e & 1)
Chris@102:       exp += restricted_multiply(result, result, a, max_bits, error);
Chris@102:    return exp;
Chris@102: }
Chris@102: 
Chris@102: inline int get_round_mode(const cpp_int& what, boost::int64_t location, boost::int64_t error)
Chris@102: {
Chris@102:    //
Chris@102:    // Can we round what at /location/, if the error in what is /error/ in
Chris@102:    // units of 0.5ulp.  Return:
Chris@102:    //
Chris@102:    // -1: Can't round.
Chris@102:    //  0: leave as is.
Chris@102:    //  1: tie.
Chris@102:    //  2: round up.
Chris@102:    //
Chris@102:    BOOST_ASSERT(location >= 0);
Chris@102:    BOOST_ASSERT(location < INT_MAX);
Chris@102:    boost::int64_t error_radius = error & 1 ? (1 + error) / 2 : error / 2;
Chris@102:    if(error_radius && ((int)msb(error_radius) >= location))
Chris@102:       return -1;
Chris@102:    if(bit_test(what, static_cast<unsigned>(location)))
Chris@102:    {
Chris@102:       if((int)lsb(what) == location)
Chris@102:          return error ? -1 : 1;   // Either a tie or can't round depending on whether we have any error
Chris@102:       if(!error)
Chris@102:          return 2;  // no error, round up.
Chris@102:       cpp_int t = what - error_radius;
Chris@102:       if((int)lsb(t) >= location)
Chris@102:          return -1;
Chris@102:       return 2;
Chris@102:    }
Chris@102:    else if(error)
Chris@102:    {
Chris@102:       cpp_int t = what + error_radius;
Chris@102:       return bit_test(t, static_cast<unsigned>(location)) ? -1 : 0;
Chris@102:    }
Chris@102:    return 0;
Chris@102: }
Chris@102: 
Chris@102: inline int get_round_mode(cpp_int& r, cpp_int& d, boost::int64_t error, const cpp_int& q)
Chris@102: {
Chris@102:    //
Chris@102:    // Lets suppose we have an inexact division by d+delta, where the true
Chris@102:    // value for the divisor is d, and with |delta| <= error/2, then
Chris@102:    // we have calculated q and r such that:
Chris@102:    //
Chris@102:    // n                  r
Chris@102:    // ---       = q + -----------
Chris@102:    // d + error        d + error
Chris@102:    //
Chris@102:    // Rearranging for n / d we get:
Chris@102:    //
Chris@102:    //    n         delta*q + r
Chris@102:    //   --- = q + -------------
Chris@102:    //    d              d
Chris@102:    //
Chris@102:    // So rounding depends on whether 2r + error * q > d.
Chris@102:    //
Chris@102:    // We return:
Chris@102:    //  0 = down down.
Chris@102:    //  1 = tie.
Chris@102:    //  2 = round up.
Chris@102:    // -1 = couldn't decide.
Chris@102:    //
Chris@102:    r <<= 1;
Chris@102:    int c = r.compare(d);
Chris@102:    if(c == 0)
Chris@102:       return error ? -1 : 1;
Chris@102:    if(c > 0)
Chris@102:    {
Chris@102:       if(error)
Chris@102:       {
Chris@102:          r -= error * q;
Chris@102:          return r.compare(d) > 0 ? 2 : -1;
Chris@102:       }
Chris@102:       return 2;
Chris@102:    }
Chris@102:    if(error)
Chris@102:    {
Chris@102:       r += error * q;
Chris@102:       return r.compare(d) < 0 ? 0 : -1;
Chris@102:    }
Chris@102:    return 0;
Chris@102: }
Chris@102: 
Chris@102: } // namespace
Chris@102: 
Chris@102: namespace backends{
Chris@102: 
Chris@102: template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
Chris@102: cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::operator=(const char *s)
Chris@102: {
Chris@102:    cpp_int n;
Chris@102:    boost::intmax_t decimal_exp = 0;
Chris@102:    boost::intmax_t digits_seen = 0;
Chris@102:    static const boost::intmax_t max_digits_seen = 4 + (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count * 301L) / 1000;
Chris@102:    bool ss = false;
Chris@102:    //
Chris@102:    // Extract the sign:
Chris@102:    //
Chris@102:    if(*s == '-')
Chris@102:    {
Chris@102:       ss = true;
Chris@102:       ++s;
Chris@102:    }
Chris@102:    else if(*s == '+')
Chris@102:       ++s;
Chris@102:    //
Chris@102:    // Special cases first:
Chris@102:    //
Chris@102:    if((std::strcmp(s, "nan") == 0) || (std::strcmp(s, "NaN") == 0) || (std::strcmp(s, "NAN") == 0))
Chris@102:    {
Chris@102:       return *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
Chris@102:    }
Chris@102:    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))
Chris@102:    {
Chris@102:       *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
Chris@102:       if(ss)
Chris@102:          negate();
Chris@102:       return *this;
Chris@102:    }
Chris@102:    //
Chris@102:    // Digits before the point:
Chris@102:    //
Chris@102:    while(*s && (*s >= '0') && (*s <= '9'))
Chris@102:    {
Chris@102:       n *= 10u;
Chris@102:       n += *s - '0';
Chris@102:       if(digits_seen || (*s != '0'))
Chris@102:          ++digits_seen;
Chris@102:       ++s;
Chris@102:    }
Chris@102:    // The decimal point (we really should localise this!!)
Chris@102:    if(*s && (*s == '.'))
Chris@102:       ++s;
Chris@102:    //
Chris@102:    // Digits after the point:
Chris@102:    //
Chris@102:    while(*s && (*s >= '0') && (*s <= '9'))
Chris@102:    {
Chris@102:       n *= 10u;
Chris@102:       n += *s - '0';
Chris@102:       --decimal_exp;
Chris@102:       if(digits_seen || (*s != '0'))
Chris@102:          ++digits_seen;
Chris@102:       ++s;
Chris@102:       if(digits_seen > max_digits_seen)
Chris@102:          break;
Chris@102:    }
Chris@102:    //
Chris@102:    // Digits we're skipping:
Chris@102:    //
Chris@102:    while(*s && (*s >= '0') && (*s <= '9'))
Chris@102:       ++s;
Chris@102:    //
Chris@102:    // See if there's an exponent:
Chris@102:    //
Chris@102:    if(*s && ((*s == 'e') || (*s == 'E')))
Chris@102:    {
Chris@102:       ++s;
Chris@102:       boost::intmax_t e = 0;
Chris@102:       bool es = false;
Chris@102:       if(*s && (*s == '-'))
Chris@102:       {
Chris@102:          es = true;
Chris@102:          ++s;
Chris@102:       }
Chris@102:       else if(*s && (*s == '+'))
Chris@102:          ++s;
Chris@102:       while(*s && (*s >= '0') && (*s <= '9'))
Chris@102:       {
Chris@102:          e *= 10u;
Chris@102:          e += *s - '0';
Chris@102:          ++s;
Chris@102:       }
Chris@102:       if(es)
Chris@102:          e = -e;
Chris@102:       decimal_exp += e;
Chris@102:    }
Chris@102:    if(*s)
Chris@102:    {
Chris@102:       //
Chris@102:       // Oops unexpected input at the end of the number:
Chris@102:       //
Chris@102:       BOOST_THROW_EXCEPTION(std::runtime_error("Unable to parse string as a valid floating point number."));
Chris@102:    }
Chris@102:    if(n == 0)
Chris@102:    {
Chris@102:       // Result is necessarily zero:
Chris@102:       *this = static_cast<limb_type>(0u);
Chris@102:       return *this;
Chris@102:    }
Chris@102: 
Chris@102:    static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
Chris@102:    //
Chris@102:    // Set our working precision - this is heuristic based, we want
Chris@102:    // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations
Chris@102:    // and excessive memory usage, but we also want to avoid having to
Chris@102:    // up the computation and start again at a higher precision.
Chris@102:    // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add
Chris@102:    // one limb for good measure.  This works very well for small exponents,
Chris@102:    // but for larger exponents we may may need to restart, we could add some
Chris@102:    // extra precision right from the start for larger exponents, but this
Chris@102:    // seems to be slightly slower in the *average* case:
Chris@102:    //
Chris@102: #ifdef BOOST_MP_STRESS_IO
Chris@102:    boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32;
Chris@102: #else
Chris@102:    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;
Chris@102: #endif
Chris@102:    boost::int64_t error = 0;
Chris@102:    boost::intmax_t calc_exp = 0;
Chris@102:    boost::intmax_t final_exponent = 0;
Chris@102: 
Chris@102:    if(decimal_exp >= 0)
Chris@102:    {
Chris@102:       // Nice and simple, the result is an integer...
Chris@102:       do
Chris@102:       {
Chris@102:          cpp_int t;
Chris@102:          if(decimal_exp)
Chris@102:          {
Chris@102:             calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), decimal_exp, max_bits, error);
Chris@102:             calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(t, t, n, max_bits, error);
Chris@102:          }
Chris@102:          else
Chris@102:             t = n;
Chris@102:          final_exponent = (boost::int64_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp + calc_exp;
Chris@102:          int rshift = msb(t) - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
Chris@102:          if(rshift > 0)
Chris@102:          {
Chris@102:             final_exponent += rshift;
Chris@102:             int roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(t, rshift - 1, error);
Chris@102:             t >>= rshift;
Chris@102:             if((roundup == 2) || ((roundup == 1) && t.backend().limbs()[0] & 1))
Chris@102:                ++t;
Chris@102:             else if(roundup < 0)
Chris@102:             {
Chris@102: #ifdef BOOST_MP_STRESS_IO
Chris@102:                max_bits += 32;
Chris@102: #else
Chris@102:                max_bits *= 2;
Chris@102: #endif
Chris@102:                error = 0;
Chris@102:                continue;
Chris@102:             }
Chris@102:          }
Chris@102:          else
Chris@102:          {
Chris@102:             BOOST_ASSERT(!error);
Chris@102:          }
Chris@102:          if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
Chris@102:          {
Chris@102:             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
Chris@102:             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
Chris@102:          }
Chris@102:          else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
Chris@102:          {
Chris@102:             // Underflow:
Chris@102:             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
Chris@102:             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
Chris@102:          }
Chris@102:          else
Chris@102:          {
Chris@102:             exponent() = static_cast<Exponent>(final_exponent);
Chris@102:             final_exponent = 0;
Chris@102:          }
Chris@102:          copy_and_round(*this, t.backend());
Chris@102:          break;
Chris@102:       }
Chris@102:       while(true);
Chris@102: 
Chris@102:       if(ss != sign())
Chris@102:          negate();
Chris@102:    }
Chris@102:    else
Chris@102:    {
Chris@102:       // Result is the ratio of two integers: we need to organise the
Chris@102:       // division so as to produce at least an N-bit result which we can
Chris@102:       // round according to the remainder.
Chris@102:       //cpp_int d = pow(cpp_int(5), -decimal_exp);
Chris@102:       do
Chris@102:       {
Chris@102:          cpp_int d;
Chris@102:          calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -decimal_exp, max_bits, error);
Chris@102:          int shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb(n) + msb(d);
Chris@102:          final_exponent = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp - calc_exp;
Chris@102:          if(shift > 0)
Chris@102:          {
Chris@102:             n <<= shift;
Chris@102:             final_exponent -= static_cast<Exponent>(shift);
Chris@102:          }
Chris@102:          cpp_int q, r;
Chris@102:          divide_qr(n, d, q, r);
Chris@102:          int gb = msb(q);
Chris@102:          BOOST_ASSERT((gb >= static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - 1));
Chris@102:          //
Chris@102:          // Check for rounding conditions we have to
Chris@102:          // handle ourselves:
Chris@102:          //
Chris@102:          int roundup = 0;
Chris@102:          if(gb == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
Chris@102:          {
Chris@102:             // Exactly the right number of bits, use the remainder to round:
Chris@102:             roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, q);
Chris@102:          }
Chris@102:          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)))
Chris@102:          {
Chris@102:             // Too many bits in q and the bits in q indicate a tie, but we can break that using r,
Chris@102:             // note that the radius of error in r is error/2 * q:
Chris@102:             int shift = gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
Chris@102:             q >>= shift;
Chris@102:             final_exponent += static_cast<Exponent>(shift);
Chris@102:             BOOST_ASSERT((msb(q) >= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
Chris@102:             if(error && (r < (error / 2) * q))
Chris@102:                roundup = -1;
Chris@102:             else if(error && (r + (error / 2) * q >= d))
Chris@102:                roundup = -1;
Chris@102:             else
Chris@102:                roundup = r ? 2 : 1;
Chris@102:          }
Chris@102:          else if(error && (((error / 2) * q + r >= d) || (r < (error / 2) * q)))
Chris@102:          {
Chris@102:             // We might have been rounding up, or got the wrong quotient: can't tell!
Chris@102:             roundup = -1;
Chris@102:          }
Chris@102:          if(roundup < 0)
Chris@102:          {
Chris@102: #ifdef BOOST_MP_STRESS_IO
Chris@102:             max_bits += 32;
Chris@102: #else
Chris@102:             max_bits *= 2;
Chris@102: #endif
Chris@102:             error = 0;
Chris@102:             if(shift > 0)
Chris@102:             {
Chris@102:                n >>= shift;
Chris@102:                final_exponent += static_cast<Exponent>(shift);
Chris@102:             }
Chris@102:             continue;
Chris@102:          }
Chris@102:          else if((roundup == 2) || ((roundup == 1) && q.backend().limbs()[0] & 1))
Chris@102:             ++q;
Chris@102:          if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
Chris@102:          {
Chris@102:             // Overflow:
Chris@102:             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
Chris@102:             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
Chris@102:          }
Chris@102:          else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
Chris@102:          {
Chris@102:             // Underflow:
Chris@102:             exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
Chris@102:             final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
Chris@102:          }
Chris@102:          else
Chris@102:          {
Chris@102:             exponent() = static_cast<Exponent>(final_exponent);
Chris@102:             final_exponent = 0;
Chris@102:          }
Chris@102:          copy_and_round(*this, q.backend());
Chris@102:          if(ss != sign())
Chris@102:             negate();
Chris@102:          break;
Chris@102:       }
Chris@102:       while(true);
Chris@102:    }
Chris@102:    //
Chris@102:    // Check for scaling and/or over/under-flow:
Chris@102:    //
Chris@102:    final_exponent += exponent();
Chris@102:    if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
Chris@102:    {
Chris@102:       // Overflow:
Chris@102:       exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
Chris@102:       bits() = limb_type(0);
Chris@102:    }
Chris@102:    else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
Chris@102:    {
Chris@102:       // Underflow:
Chris@102:       exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
Chris@102:       bits() = limb_type(0);
Chris@102:       sign() = 0;
Chris@102:    }
Chris@102:    else
Chris@102:    {
Chris@102:       exponent() = static_cast<Exponent>(final_exponent);
Chris@102:    }
Chris@102:    return *this;
Chris@102: }
Chris@102: 
Chris@102: template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
Chris@102: std::string cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::str(std::streamsize dig, std::ios_base::fmtflags f) const
Chris@102: {
Chris@102:    if(dig == 0)
Chris@102:       dig = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::max_digits10;
Chris@102: 
Chris@102:    bool scientific = (f & std::ios_base::scientific) == std::ios_base::scientific;
Chris@102:    bool fixed = !scientific && (f & std::ios_base::fixed);
Chris@102: 
Chris@102:    std::string s;
Chris@102: 
Chris@102:    if(exponent() <= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
Chris@102:    {
Chris@102:       // How far to left-shift in order to demormalise the mantissa:
Chris@102:       boost::intmax_t shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1;
Chris@102:       boost::intmax_t digits_wanted = static_cast<int>(dig);
Chris@102:       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()));
Chris@102:       //
Chris@102:       // For fixed formatting we want /dig/ digits after the decimal point,
Chris@102:       // so if the exponent is zero, allowing for the one digit before the
Chris@102:       // decimal point, we want 1 + dig digits etc.
Chris@102:       //
Chris@102:       if(fixed)
Chris@102:          digits_wanted += 1 + base10_exp;
Chris@102:       if(scientific)
Chris@102:          digits_wanted += 1;
Chris@102:       if(digits_wanted < -1)
Chris@102:       {
Chris@102:          // Fixed precision, no significant digits, and nothing to round!
Chris@102:          s = "0";
Chris@102:          if(sign())
Chris@102:             s.insert(static_cast<std::string::size_type>(0), 1, '-');
Chris@102:          boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, true);
Chris@102:          return s;
Chris@102:       }
Chris@102:       //
Chris@102:       // power10 is the base10 exponent we need to multiply/divide by in order
Chris@102:       // to convert our denormalised number to an integer with the right number of digits:
Chris@102:       //
Chris@102:       boost::intmax_t power10 = digits_wanted - base10_exp - 1;
Chris@102:       //
Chris@102:       // If we calculate 5^power10 rather than 10^power10 we need to move
Chris@102:       // 2^power10 into /shift/
Chris@102:       //
Chris@102:       shift -= power10;
Chris@102:       cpp_int i;
Chris@102:       int roundup = 0; // 0=no rounding, 1=tie, 2=up
Chris@102:       static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
Chris@102:       //
Chris@102:       // Set our working precision - this is heuristic based, we want
Chris@102:       // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations
Chris@102:       // and excessive memory usage, but we also want to avoid having to
Chris@102:       // up the computation and start again at a higher precision.
Chris@102:       // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add
Chris@102:       // one limb for good measure.  This works very well for small exponents,
Chris@102:       // but for larger exponents we add a few extra limbs to max_bits:
Chris@102:       //
Chris@102: #ifdef BOOST_MP_STRESS_IO
Chris@102:       boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32;
Chris@102: #else
Chris@102:       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;
Chris@102:       if(power10)
Chris@102:          max_bits += (msb(boost::multiprecision::detail::abs(power10)) / 8) * limb_bits;
Chris@102: #endif
Chris@102:       do
Chris@102:       {
Chris@102:          boost::int64_t error = 0;
Chris@102:          boost::intmax_t calc_exp = 0;
Chris@102:          //
Chris@102:          // Our integer result is: bits() * 2^-shift * 5^power10
Chris@102:          //
Chris@102:          i = bits();
Chris@102:          if(shift < 0)
Chris@102:          {
Chris@102:             if(power10 >= 0)
Chris@102:             {
Chris@102:                // We go straight to the answer with all integer arithmetic,
Chris@102:                // the result is always exact and never needs rounding:
Chris@102:                BOOST_ASSERT(power10 <= (boost::intmax_t)INT_MAX);
Chris@102:                i <<= -shift;
Chris@102:                if(power10)
Chris@102:                   i *= pow(cpp_int(5), static_cast<unsigned>(power10));
Chris@102:             }
Chris@102:             else if(power10 < 0)
Chris@102:             {
Chris@102:                cpp_int d;
Chris@102:                calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -power10, max_bits, error);
Chris@102:                shift += calc_exp;
Chris@102:                BOOST_ASSERT(shift < 0); // Must still be true!
Chris@102:                i <<= -shift;
Chris@102:                cpp_int r;
Chris@102:                divide_qr(i, d, i, r);
Chris@102:                roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, i);
Chris@102:                if(roundup < 0)
Chris@102:                {
Chris@102: #ifdef BOOST_MP_STRESS_IO
Chris@102:                   max_bits += 32;
Chris@102: #else
Chris@102:                   max_bits *= 2;
Chris@102: #endif
Chris@102:                   shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
Chris@102:                   continue;
Chris@102:                }
Chris@102:             }
Chris@102:          }
Chris@102:          else
Chris@102:          {
Chris@102:             //
Chris@102:             // Our integer is bits() * 2^-shift * 10^power10
Chris@102:             //
Chris@102:             if(power10 > 0)
Chris@102:             {
Chris@102:                if(power10)
Chris@102:                {
Chris@102:                   cpp_int t;
Chris@102:                   calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), power10, max_bits, error);
Chris@102:                   calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(i, i, t, max_bits, error);
Chris@102:                   shift -= calc_exp;
Chris@102:                }
Chris@102:                if((shift < 0) || ((shift == 0) && error))
Chris@102:                {
Chris@102:                   // We only get here if we were asked for a crazy number of decimal digits -
Chris@102:                   // more than are present in a 2^max_bits number.
Chris@102: #ifdef BOOST_MP_STRESS_IO
Chris@102:                   max_bits += 32;
Chris@102: #else
Chris@102:                   max_bits *= 2;
Chris@102: #endif
Chris@102:                   shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
Chris@102:                   continue;
Chris@102:                }
Chris@102:                if(shift)
Chris@102:                {
Chris@102:                   roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(i, shift - 1, error);
Chris@102:                   if(roundup < 0)
Chris@102:                   {
Chris@102: #ifdef BOOST_MP_STRESS_IO
Chris@102:                      max_bits += 32;
Chris@102: #else
Chris@102:                      max_bits *= 2;
Chris@102: #endif
Chris@102:                      shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
Chris@102:                      continue;
Chris@102:                   }
Chris@102:                   i >>= shift;
Chris@102:                }
Chris@102:             }
Chris@102:             else
Chris@102:             {
Chris@102:                // We're right shifting, *and* dividing by 5^-power10,
Chris@102:                // so 5^-power10 can never be that large or we'd simply
Chris@102:                // get zero as a result, and that case is already handled above:
Chris@102:                cpp_int r;
Chris@102:                BOOST_ASSERT(-power10 < INT_MAX);
Chris@102:                cpp_int d = pow(cpp_int(5), static_cast<unsigned>(-power10));
Chris@102:                d <<= shift;
Chris@102:                divide_qr(i, d, i, r);
Chris@102:                r <<= 1;
Chris@102:                int c = r.compare(d);
Chris@102:                roundup = c < 0 ? 0 : c == 0 ? 1 : 2;
Chris@102:             }
Chris@102:          }
Chris@102:          s = i.str(0, std::ios_base::fmtflags(0));
Chris@102:          //
Chris@102:          // Check if we got the right number of digits, this
Chris@102:          // is really a test of whether we calculated the
Chris@102:          // decimal exponent correctly:
Chris@102:          //
Chris@102:          boost::intmax_t digits_got = i ? static_cast<boost::intmax_t>(s.size()) : 0;
Chris@102:          if(digits_got != digits_wanted)
Chris@102:          {
Chris@102:             base10_exp += digits_got - digits_wanted;
Chris@102:             if(fixed)
Chris@102:                digits_wanted = digits_got;  // strange but true.
Chris@102:             power10 = digits_wanted - base10_exp - 1;
Chris@102:             shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
Chris@102:             if(fixed)
Chris@102:                break;
Chris@102:             roundup = 0;
Chris@102:          }
Chris@102:          else
Chris@102:             break;
Chris@102:       }
Chris@102:       while(true);
Chris@102:       //
Chris@102:       // Check whether we need to round up: note that we could equally round up
Chris@102:       // the integer /i/ above, but since we need to perform the rounding *after*
Chris@102:       // the conversion to a string and the digit count check, we might as well
Chris@102:       // do it here:
Chris@102:       //
Chris@102:       if((roundup == 2) || ((roundup == 1) && ((s[s.size() - 1] - '0') & 1)))
Chris@102:       {
Chris@102:          boost::multiprecision::detail::round_string_up_at(s, static_cast<int>(s.size() - 1), base10_exp);
Chris@102:       }
Chris@102: 
Chris@102:       if(sign())
Chris@102:          s.insert(static_cast<std::string::size_type>(0), 1, '-');
Chris@102: 
Chris@102:       boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, false);
Chris@102:    }
Chris@102:    else
Chris@102:    {
Chris@102:       switch(exponent())
Chris@102:       {
Chris@102:       case exponent_zero:
Chris@102:          s = "0";
Chris@102:          boost::multiprecision::detail::format_float_string(s, 0, dig, f, true);
Chris@102:          break;
Chris@102:       case exponent_nan:
Chris@102:          s = "nan";
Chris@102:          break;
Chris@102:       case exponent_infinity:
Chris@102:          s = sign() ? "-inf" : f & std::ios_base::showpos ? "+inf" : "inf";
Chris@102:          break;
Chris@102:       }
Chris@102:    }
Chris@102:    return s;
Chris@102: }
Chris@102: 
Chris@102: }}} // namespaces
Chris@102: 
Chris@102: #endif
Chris@102: