boost/multiprecision/cpp_bin_float.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_MATH_CPP_BIN_FLOAT_HPP
   7 #define BOOST_MATH_CPP_BIN_FLOAT_HPP
   8
   9 #include <boost/multiprecision/cpp_int.hpp>
  10 #include <boost/multiprecision/integer.hpp>
  11 #include <boost/math/special_functions/trunc.hpp>
  12 #include <boost/multiprecision/detail/float_string_cvt.hpp>
  13
  14 //
  15 // Some includes we need from Boost.Math, since we rely on that library to provide these functions:
  16 //
  17 #include <boost/math/special_functions/asinh.hpp>
  18 #include <boost/math/special_functions/acosh.hpp>
  19 #include <boost/math/special_functions/atanh.hpp>
  20 #include <boost/math/special_functions/cbrt.hpp>
  21 #include <boost/math/special_functions/expm1.hpp>
  22 #include <boost/math/special_functions/gamma.hpp>
  23
  24 #ifdef BOOST_HAS_FLOAT128
  25 #include <quadmath.h>
  26 #endif
  27
  28 namespace boost {
  29 namespace multiprecision {
  30 namespace backends {
  31
  32 enum digit_base_type
  33 {
  34    digit_base_2  = 2,
  35    digit_base_10 = 10
  36 };
  37
  38 #ifdef BOOST_MSVC
  39 #pragma warning(push)
  40 #pragma warning(disable : 4522 6326) // multiple assignment operators specified, comparison of two constants
  41 #endif
  42
  43 namespace detail {
  44
  45 template <class U>
  46 inline typename enable_if_c<is_unsigned<U>::value, bool>::type is_negative(U) { return false; }
  47 template <class S>
  48 inline typename disable_if_c<is_unsigned<S>::value, bool>::type is_negative(S s) { return s < 0; }
  49
  50 template <class Float, int, bool = number_category<Float>::value == number_kind_floating_point>
  51 struct is_cpp_bin_float_implicitly_constructible_from_type
  52 {
  53    static const bool value = false;
  54 };
  55
  56 template <class Float, int bit_count>
  57 struct is_cpp_bin_float_implicitly_constructible_from_type<Float, bit_count, true>
  58 {
  59    static const bool value = (std::numeric_limits<Float>::digits <= (int)bit_count) && (std::numeric_limits<Float>::radix == 2) && std::numeric_limits<Float>::is_specialized
  60 #ifdef BOOST_HAS_FLOAT128
  61                              && !boost::is_same<Float, __float128>::value
  62 #endif
  63                              && (is_floating_point<Float>::value || is_number<Float>::value);
  64 };
  65
  66 template <class Float, int, bool = number_category<Float>::value == number_kind_floating_point>
  67 struct is_cpp_bin_float_explicitly_constructible_from_type
  68 {
  69    static const bool value = false;
  70 };
  71
  72 template <class Float, int bit_count>
  73 struct is_cpp_bin_float_explicitly_constructible_from_type<Float, bit_count, true>
  74 {
  75    static const bool value = (std::numeric_limits<Float>::digits > (int)bit_count) && (std::numeric_limits<Float>::radix == 2) && std::numeric_limits<Float>::is_specialized
  76 #ifdef BOOST_HAS_FLOAT128
  77                              && !boost::is_same<Float, __float128>::value
  78 #endif
  79        ;
  80 };
  81
  82 } // namespace detail
  83
  84 template <unsigned Digits, digit_base_type DigitBase = digit_base_10, class Allocator = void, class Exponent = int, Exponent MinExponent = 0, Exponent MaxExponent = 0>
  85 class cpp_bin_float
  86 {
  87  public:
  88    static const unsigned                                                                                                                                                          bit_count = DigitBase == digit_base_2 ? Digits : (Digits * 1000uL) / 301uL + (((Digits * 1000uL) % 301) ? 2u : 1u);
  89    typedef cpp_int_backend<is_void<Allocator>::value ? bit_count : 0, bit_count, is_void<Allocator>::value ? unsigned_magnitude : signed_magnitude, unchecked, Allocator>         rep_type;
  90    typedef cpp_int_backend<is_void<Allocator>::value ? 2 * bit_count : 0, 2 * bit_count, is_void<Allocator>::value ? unsigned_magnitude : signed_magnitude, unchecked, Allocator> double_rep_type;
  91
  92    typedef typename rep_type::signed_types              signed_types;
  93    typedef typename rep_type::unsigned_types            unsigned_types;
  94    typedef boost::mpl::list<float, double, long double> float_types;
  95    typedef Exponent                                     exponent_type;
  96
  97    static const exponent_type max_exponent_limit = boost::integer_traits<exponent_type>::const_max - 2 * static_cast<exponent_type>(bit_count);
  98    static const exponent_type min_exponent_limit = boost::integer_traits<exponent_type>::const_min + 2 * static_cast<exponent_type>(bit_count);
  99
 100    BOOST_STATIC_ASSERT_MSG(MinExponent >= min_exponent_limit, "Template parameter MinExponent is too negative for our internal logic to function correctly, sorry!");
 101    BOOST_STATIC_ASSERT_MSG(MaxExponent <= max_exponent_limit, "Template parameter MaxExponent is too large for our internal logic to function correctly, sorry!");
 102    BOOST_STATIC_ASSERT_MSG(MinExponent <= 0, "Template parameter MinExponent can not be positive!");
 103    BOOST_STATIC_ASSERT_MSG(MaxExponent >= 0, "Template parameter MaxExponent can not be negative!");
 104
 105    static const exponent_type max_exponent = MaxExponent == 0 ? max_exponent_limit : MaxExponent;
 106    static const exponent_type min_exponent = MinExponent == 0 ? min_exponent_limit : MinExponent;
 107
 108    static const exponent_type exponent_zero     = max_exponent + 1;
 109    static const exponent_type exponent_infinity = max_exponent + 2;
 110    static const exponent_type exponent_nan      = max_exponent + 3;
 111
 112  private:
 113    rep_type      m_data;
 114    exponent_type m_exponent;
 115    bool          m_sign;
 116
 117  public:
 118    cpp_bin_float() BOOST_MP_NOEXCEPT_IF(noexcept(rep_type())) : m_data(), m_exponent(exponent_zero), m_sign(false) {}
 119
 120    cpp_bin_float(const cpp_bin_float& o) BOOST_MP_NOEXCEPT_IF(noexcept(rep_type(std::declval<const rep_type&>())))
 121        : m_data(o.m_data), m_exponent(o.m_exponent), m_sign(o.m_sign) {}
 122
 123    template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
 124    cpp_bin_float(const cpp_bin_float<D, B, A, E, MinE, MaxE>& o, typename boost::enable_if_c<(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = 0)
 125    {
 126       *this = o;
 127    }
 128    template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
 129    explicit cpp_bin_float(const cpp_bin_float<D, B, A, E, MinE, MaxE>& o, typename boost::disable_if_c<(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = 0)
 130        : m_exponent(o.exponent()), m_sign(o.sign())
 131    {
 132       *this = o;
 133    }
 134    template <class Float>
 135    cpp_bin_float(const Float& f,
 136                  typename boost::enable_if_c<detail::is_cpp_bin_float_implicitly_constructible_from_type<Float, bit_count>::value>::type const* = 0)
 137        : m_data(), m_exponent(0), m_sign(false)
 138    {
 139       this->assign_float(f);
 140    }
 141
 142    template <class Float>
 143    explicit cpp_bin_float(const Float& f,
 144                           typename boost::enable_if_c<detail::is_cpp_bin_float_explicitly_constructible_from_type<Float, bit_count>::value>::type const* = 0)
 145        : m_data(), m_exponent(0), m_sign(false)
 146    {
 147       this->assign_float(f);
 148    }
 149 #ifdef BOOST_HAS_FLOAT128
 150    template <class Float>
 151    cpp_bin_float(const Float& f,
 152                  typename boost::enable_if_c<
 153                      boost::is_same<Float, __float128>::value && ((int)bit_count >= 113)>::type const* = 0)
 154        : m_data(), m_exponent(0), m_sign(false)
 155    {
 156       this->assign_float(f);
 157    }
 158    template <class Float>
 159    explicit cpp_bin_float(const Float& f,
 160                           typename boost::enable_if_c<
 161                               boost::is_same<Float, __float128>::value && ((int)bit_count < 113)>::type const* = 0)
 162        : m_data(), m_exponent(0), m_sign(false)
 163    {
 164       this->assign_float(f);
 165    }
 166 #endif
 167    cpp_bin_float& operator=(const cpp_bin_float& o) BOOST_MP_NOEXCEPT_IF(noexcept(std::declval<rep_type&>() = std::declval<const rep_type&>()))
 168    {
 169       m_data     = o.m_data;
 170       m_exponent = o.m_exponent;
 171       m_sign     = o.m_sign;
 172       return *this;
 173    }
 174
 175    template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
 176    cpp_bin_float& operator=(const cpp_bin_float<D, B, A, E, MinE, MaxE>& f)
 177    {
 178       switch (eval_fpclassify(f))
 179       {
 180       case FP_ZERO:
 181          m_data     = limb_type(0);
 182          m_sign     = f.sign();
 183          m_exponent = exponent_zero;
 184          break;
 185       case FP_NAN:
 186          m_data     = limb_type(0);
 187          m_sign     = false;
 188          m_exponent = exponent_nan;
 189          break;
 190          ;
 191       case FP_INFINITE:
 192          m_data     = limb_type(0);
 193          m_sign     = f.sign();
 194          m_exponent = exponent_infinity;
 195          break;
 196       default:
 197          typename cpp_bin_float<D, B, A, E, MinE, MaxE>::rep_type b(f.bits());
 198          this->exponent() = f.exponent() + (E)bit_count - (E)cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count;
 199          this->sign()     = f.sign();
 200          copy_and_round(*this, b);
 201       }
 202       return *this;
 203    }
 204 #ifdef BOOST_HAS_FLOAT128
 205    template <class Float>
 206    typename boost::enable_if_c<
 207        (number_category<Float>::value == number_kind_floating_point)
 208            //&& (std::numeric_limits<Float>::digits <= (int)bit_count)
 209            && ((std::numeric_limits<Float>::radix == 2) || (boost::is_same<Float, __float128>::value)),
 210        cpp_bin_float&>::type
 211    operator=(const Float& f)
 212 #else
 213    template <class Float>
 214    typename boost::enable_if_c<
 215        (number_category<Float>::value == number_kind_floating_point)
 216            //&& (std::numeric_limits<Float>::digits <= (int)bit_count)
 217            && (std::numeric_limits<Float>::radix == 2),
 218        cpp_bin_float&>::type
 219    operator=(const Float& f)
 220 #endif
 221    {
 222       return assign_float(f);
 223    }
 224
 225 #ifdef BOOST_HAS_FLOAT128
 226    template <class Float>
 227    typename boost::enable_if_c<boost::is_same<Float, __float128>::value, cpp_bin_float&>::type assign_float(Float f)
 228    {
 229       using default_ops::eval_add;
 230       typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type;
 231       if (f == 0)
 232       {
 233          m_data     = limb_type(0);
 234          m_sign     = (signbitq(f) > 0);
 235          m_exponent = exponent_zero;
 236          return *this;
 237       }
 238       else if (isnanq(f))
 239       {
 240          m_data     = limb_type(0);
 241          m_sign     = false;
 242          m_exponent = exponent_nan;
 243          return *this;
 244       }
 245       else if (isinfq(f))
 246       {
 247          m_data     = limb_type(0);
 248          m_sign     = (f < 0);
 249          m_exponent = exponent_infinity;
 250          return *this;
 251       }
 252       if (f < 0)
 253       {
 254          *this = -f;
 255          this->negate();
 256          return *this;
 257       }
 258
 259       typedef typename mpl::front<unsigned_types>::type ui_type;
 260       m_data     = static_cast<ui_type>(0u);
 261       m_sign     = false;
 262       m_exponent = 0;
 263
 264       static const int bits = sizeof(int) * CHAR_BIT - 1;
 265       int              e;
 266       f = frexpq(f, &e);
 267       while (f)
 268       {
 269          f = ldexpq(f, bits);
 270          e -= bits;
 271          int ipart = (int)truncq(f);
 272          f -= ipart;
 273          m_exponent += bits;
 274          cpp_bin_float t;
 275          t = static_cast<bf_int_type>(ipart);
 276          eval_add(*this, t);
 277       }
 278       m_exponent += static_cast<Exponent>(e);
 279       return *this;
 280    }
 281 #endif
 282 #ifdef BOOST_HAS_FLOAT128
 283    template <class Float>
 284    typename boost::enable_if_c<is_floating_point<Float>::value && !is_same<Float, __float128>::value, cpp_bin_float&>::type assign_float(Float f)
 285 #else
 286    template <class Float>
 287    typename boost::enable_if_c<is_floating_point<Float>::value, cpp_bin_float&>::type assign_float(Float f)
 288 #endif
 289    {
 290       BOOST_MATH_STD_USING
 291       using default_ops::eval_add;
 292       typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type;
 293
 294       switch ((boost::math::fpclassify)(f))
 295       {
 296       case FP_ZERO:
 297          m_data     = limb_type(0);
 298          m_sign     = ((boost::math::signbit)(f) > 0);
 299          m_exponent = exponent_zero;
 300          return *this;
 301       case FP_NAN:
 302          m_data     = limb_type(0);
 303          m_sign     = false;
 304          m_exponent = exponent_nan;
 305          return *this;
 306       case FP_INFINITE:
 307          m_data     = limb_type(0);
 308          m_sign     = (f < 0);
 309          m_exponent = exponent_infinity;
 310          return *this;
 311       }
 312       if (f < 0)
 313       {
 314          *this = -f;
 315          this->negate();
 316          return *this;
 317       }
 318
 319       typedef typename mpl::front<unsigned_types>::type ui_type;
 320       m_data     = static_cast<ui_type>(0u);
 321       m_sign     = false;
 322       m_exponent = 0;
 323
 324       static const int bits = sizeof(int) * CHAR_BIT - 1;
 325       int              e;
 326       f = frexp(f, &e);
 327       while (f)
 328       {
 329          f = ldexp(f, bits);
 330          e -= bits;
 331 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 332          int ipart = itrunc(f);
 333 #else
 334          int ipart = static_cast<int>(f);
 335 #endif
 336          f -= ipart;
 337          m_exponent += bits;
 338          cpp_bin_float t;
 339          t = static_cast<bf_int_type>(ipart);
 340          eval_add(*this, t);
 341       }
 342       m_exponent += static_cast<Exponent>(e);
 343       return *this;
 344    }
 345
 346    template <class Float>
 347    typename boost::enable_if_c<
 348        (number_category<Float>::value == number_kind_floating_point) && !boost::is_floating_point<Float>::value && (number_category<Float>::value == number_kind_floating_point),
 349        cpp_bin_float&>::type
 350    assign_float(Float f)
 351    {
 352       BOOST_MATH_STD_USING
 353       using default_ops::eval_add;
 354       using default_ops::eval_convert_to;
 355       using default_ops::eval_get_sign;
 356       using default_ops::eval_subtract;
 357
 358       typedef typename boost::multiprecision::detail::canonical<int, Float>::type         f_int_type;
 359       typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type;
 360
 361       switch (eval_fpclassify(f))
 362       {
 363       case FP_ZERO:
 364          m_data     = limb_type(0);
 365          m_sign     = (eval_get_sign(f) > 0);
 366          m_exponent = exponent_zero;
 367          return *this;
 368       case FP_NAN:
 369          m_data     = limb_type(0);
 370          m_sign     = false;
 371          m_exponent = exponent_nan;
 372          return *this;
 373       case FP_INFINITE:
 374          m_data     = limb_type(0);
 375          m_sign     = eval_get_sign(f) < 0;
 376          m_exponent = exponent_infinity;
 377          return *this;
 378       }
 379       if (eval_get_sign(f) < 0)
 380       {
 381          f.negate();
 382          *this = f;
 383          this->negate();
 384          return *this;
 385       }
 386
 387       typedef typename mpl::front<unsigned_types>::type ui_type;
 388       m_data     = static_cast<ui_type>(0u);
 389       m_sign     = false;
 390       m_exponent = 0;
 391
 392       static const int bits = sizeof(int) * CHAR_BIT - 1;
 393       int              e;
 394       eval_frexp(f, f, &e);
 395       while (eval_get_sign(f) != 0)
 396       {
 397          eval_ldexp(f, f, bits);
 398          e -= bits;
 399          int ipart;
 400          eval_convert_to(&ipart, f);
 401          eval_subtract(f, static_cast<f_int_type>(ipart));
 402          m_exponent += bits;
 403          eval_add(*this, static_cast<bf_int_type>(ipart));
 404       }
 405       m_exponent += e;
 406       if (m_exponent > max_exponent)
 407          m_exponent = exponent_infinity;
 408       if (m_exponent < min_exponent)
 409       {
 410          m_data     = limb_type(0u);
 411          m_exponent = exponent_zero;
 412          m_sign     = (eval_get_sign(f) > 0);
 413       }
 414       else if (eval_get_sign(m_data) == 0)
 415       {
 416          m_exponent = exponent_zero;
 417          m_sign     = (eval_get_sign(f) > 0);
 418       }
 419       return *this;
 420    }
 421    template <class B, expression_template_option et>
 422    cpp_bin_float& assign_float(const number<B, et>& f)
 423    {
 424       return assign_float(f.backend());
 425    }
 426
 427    template <class I>
 428    typename boost::enable_if<is_integral<I>, cpp_bin_float&>::type operator=(const I& i)
 429    {
 430       using default_ops::eval_bit_test;
 431       if (!i)
 432       {
 433          m_data     = static_cast<limb_type>(0);
 434          m_exponent = exponent_zero;
 435          m_sign     = false;
 436       }
 437       else
 438       {
 439          typedef typename make_unsigned<I>::type                                            ui_type;
 440          ui_type                                                                            fi = static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(i));
 441          typedef typename boost::multiprecision::detail::canonical<ui_type, rep_type>::type ar_type;
 442          m_data         = static_cast<ar_type>(fi);
 443          unsigned shift = msb(fi);
 444          if (shift >= bit_count)
 445          {
 446             m_exponent = static_cast<Exponent>(shift);
 447             m_data     = static_cast<ar_type>(fi >> (shift + 1 - bit_count));
 448          }
 449          else
 450          {
 451             m_exponent = static_cast<Exponent>(shift);
 452             eval_left_shift(m_data, bit_count - shift - 1);
 453          }
 454          BOOST_ASSERT(eval_bit_test(m_data, bit_count - 1));
 455          m_sign = detail::is_negative(i);
 456       }
 457       return *this;
 458    }
 459
 460    cpp_bin_float& operator=(const char* s);
 461
 462    void swap(cpp_bin_float& o) BOOST_NOEXCEPT
 463    {
 464       m_data.swap(o.m_data);
 465       std::swap(m_exponent, o.m_exponent);
 466       std::swap(m_sign, o.m_sign);
 467    }
 468
 469    std::string str(std::streamsize dig, std::ios_base::fmtflags f) const;
 470
 471    void negate()
 472    {
 473       if (m_exponent != exponent_nan)
 474          m_sign = !m_sign;
 475    }
 476
 477    int compare(const cpp_bin_float& o) const BOOST_NOEXCEPT
 478    {
 479       if (m_sign != o.m_sign)
 480          return (m_exponent == exponent_zero) && (m_exponent == o.m_exponent) ? 0 : m_sign ? -1 : 1;
 481       int result;
 482       if (m_exponent == exponent_nan)
 483          return -1;
 484       else if (m_exponent != o.m_exponent)
 485       {
 486          if (m_exponent == exponent_zero)
 487             result = -1;
 488          else if (o.m_exponent == exponent_zero)
 489             result = 1;
 490          else
 491             result = m_exponent > o.m_exponent ? 1 : -1;
 492       }
 493       else
 494          result = m_data.compare(o.m_data);
 495       if (m_sign)
 496          result = -result;
 497       return result;
 498    }
 499    template <class A>
 500    int compare(const A& o) const BOOST_NOEXCEPT
 501    {
 502       cpp_bin_float b;
 503       b = o;
 504       return compare(b);
 505    }
 506
 507    rep_type&            bits() { return m_data; }
 508    const rep_type&      bits() const { return m_data; }
 509    exponent_type&       exponent() { return m_exponent; }
 510    const exponent_type& exponent() const { return m_exponent; }
 511    bool&                sign() { return m_sign; }
 512    const bool&          sign() const { return m_sign; }
 513    void                 check_invariants()
 514    {
 515       using default_ops::eval_bit_test;
 516       using default_ops::eval_is_zero;
 517       if ((m_exponent <= max_exponent) && (m_exponent >= min_exponent))
 518       {
 519          BOOST_ASSERT(eval_bit_test(m_data, bit_count - 1));
 520       }
 521       else
 522       {
 523          BOOST_ASSERT(m_exponent > max_exponent);
 524          BOOST_ASSERT(m_exponent <= exponent_nan);
 525          BOOST_ASSERT(eval_is_zero(m_data));
 526       }
 527    }
 528    template <class Archive>
 529    void serialize(Archive& ar, const unsigned int /*version*/)
 530    {
 531       ar& boost::make_nvp("data", m_data);
 532       ar& boost::make_nvp("exponent", m_exponent);
 533       ar& boost::make_nvp("sign", m_sign);
 534    }
 535 };
 536
 537 #ifdef BOOST_MSVC
 538 #pragma warning(pop)
 539 #endif
 540
 541 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class Int>
 542 inline void copy_and_round(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, Int& arg, int bits_to_keep = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
 543 {
 544    // Precondition: exponent of res must have been set before this function is called
 545    // as we may need to adjust it based on how many bits_to_keep in arg are set.
 546    using default_ops::eval_bit_test;
 547    using default_ops::eval_get_sign;
 548    using default_ops::eval_increment;
 549    using default_ops::eval_left_shift;
 550    using default_ops::eval_lsb;
 551    using default_ops::eval_msb;
 552    using default_ops::eval_right_shift;
 553
 554    // cancellation may have resulted in arg being all zeros:
 555    if (eval_get_sign(arg) == 0)
 556    {
 557       res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
 558       res.sign()     = false;
 559       res.bits()     = static_cast<limb_type>(0u);
 560       return;
 561    }
 562    int msb = eval_msb(arg);
 563    if (static_cast<int>(bits_to_keep) > msb + 1)
 564    {
 565       // Must have had cancellation in subtraction,
 566       // or be converting from a narrower type, so shift left:
 567       res.bits() = arg;
 568       eval_left_shift(res.bits(), bits_to_keep - msb - 1);
 569       res.exponent() -= static_cast<Exponent>(bits_to_keep - msb - 1);
 570    }
 571    else if (static_cast<int>(bits_to_keep) < msb + 1)
 572    {
 573       // We have more bits_to_keep than we need, so round as required,
 574       // first get the rounding bit:
 575       bool roundup = eval_bit_test(arg, msb - bits_to_keep);
 576       // Then check for a tie:
 577       if (roundup && (msb - bits_to_keep == (int)eval_lsb(arg)))
 578       {
 579          // Ties round towards even:
 580          if (!eval_bit_test(arg, msb - bits_to_keep + 1))
 581             roundup = false;
 582       }
 583       // Shift off the bits_to_keep we don't need:
 584       eval_right_shift(arg, msb - bits_to_keep + 1);
 585       res.exponent() += static_cast<Exponent>(msb - bits_to_keep + 1);
 586       if (roundup)
 587       {
 588          eval_increment(arg);
 589          if (bits_to_keep)
 590          {
 591             if (eval_bit_test(arg, bits_to_keep))
 592             {
 593                // This happens very very rairly, all the bits left after
 594                // truncation must be 1's and we're rounding up an order of magnitude:
 595                eval_right_shift(arg, 1u);
 596                ++res.exponent();
 597             }
 598          }
 599          else
 600          {
 601             // We get here when bits_to_keep is zero but we're rounding up,
 602             // as a result we end up with a single digit that is a 1:
 603             ++bits_to_keep;
 604          }
 605       }
 606       if (bits_to_keep != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
 607       {
 608          // Normalize result when we're rounding to fewer bits than we can hold, only happens in conversions
 609          // to narrower types:
 610          eval_left_shift(arg, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - bits_to_keep);
 611          res.exponent() -= static_cast<Exponent>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - bits_to_keep);
 612       }
 613       res.bits() = arg;
 614    }
 615    else
 616    {
 617       res.bits() = arg;
 618    }
 619    if (!bits_to_keep && !res.bits().limbs()[0])
 620    {
 621       // We're keeping zero bits and did not round up, so result is zero:
 622       res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
 623       return;
 624    }
 625    // Result must be normalized:
 626    BOOST_ASSERT(((int)eval_msb(res.bits()) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
 627
 628    if (res.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
 629    {
 630       // Overflow:
 631       res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
 632       res.bits()     = static_cast<limb_type>(0u);
 633    }
 634    else if (res.exponent() < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
 635    {
 636       // Underflow:
 637       res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
 638       res.bits()     = static_cast<limb_type>(0u);
 639    }
 640 }
 641
 642 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 643 inline void do_eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b)
 644 {
 645    if (a.exponent() < b.exponent())
 646    {
 647       bool s = a.sign();
 648       do_eval_add(res, b, a);
 649       if (res.sign() != s)
 650          res.negate();
 651       return;
 652    }
 653
 654    using default_ops::eval_add;
 655    using default_ops::eval_bit_test;
 656
 657    typedef typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type exponent_type;
 658
 659    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
 660
 661    // Special cases first:
 662    switch (a.exponent())
 663    {
 664    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
 665    {
 666       bool s     = a.sign();
 667       res        = b;
 668       res.sign() = s;
 669       return;
 670    }
 671    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
 672       if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)
 673          res = b;
 674       else
 675          res = a;
 676       return; // result is still infinite.
 677    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
 678       res = a;
 679       return; // result is still a NaN.
 680    }
 681    switch (b.exponent())
 682    {
 683    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
 684       res = a;
 685       return;
 686    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
 687       res = b;
 688       if (res.sign())
 689          res.negate();
 690       return; // result is infinite.
 691    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
 692       res = b;
 693       return; // result is a NaN.
 694    }
 695
 696    BOOST_STATIC_ASSERT(boost::integer_traits<exponent_type>::const_max - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent);
 697
 698    bool s = a.sign();
 699    dt     = a.bits();
 700    if (a.exponent() > (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent())
 701    {
 702       res.exponent() = a.exponent();
 703    }
 704    else
 705    {
 706       exponent_type e_diff = a.exponent() - b.exponent();
 707       BOOST_ASSERT(e_diff >= 0);
 708       eval_left_shift(dt, e_diff);
 709       res.exponent() = a.exponent() - e_diff;
 710       eval_add(dt, b.bits());
 711    }
 712
 713    copy_and_round(res, dt);
 714    res.check_invariants();
 715    if (res.sign() != s)
 716       res.negate();
 717 }
 718
 719 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 720 inline void do_eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b)
 721 {
 722    using default_ops::eval_bit_test;
 723    using default_ops::eval_decrement;
 724    using default_ops::eval_subtract;
 725
 726    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
 727
 728    // Special cases first:
 729    switch (a.exponent())
 730    {
 731    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
 732       if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)
 733          res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
 734       else
 735       {
 736          bool s = a.sign();
 737          res    = b;
 738          if (res.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero)
 739             res.sign() = false;
 740          else if (res.sign() == s)
 741             res.negate();
 742       }
 743       return;
 744    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
 745       if ((b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan) || (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity))
 746          res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
 747       else
 748          res = a;
 749       return;
 750    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
 751       res = a;
 752       return; // result is still a NaN.
 753    }
 754    switch (b.exponent())
 755    {
 756    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
 757       res = a;
 758       return;
 759    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
 760       res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
 761       res.sign()     = !a.sign();
 762       res.bits()     = static_cast<limb_type>(0u);
 763       return; // result is a NaN.
 764    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
 765       res = b;
 766       return; // result is still a NaN.
 767    }
 768
 769    bool s = a.sign();
 770    if ((a.exponent() > b.exponent()) || ((a.exponent() == b.exponent()) && a.bits().compare(b.bits()) >= 0))
 771    {
 772       dt = a.bits();
 773       if (a.exponent() <= (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent())
 774       {
 775          typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e_diff = a.exponent() - b.exponent();
 776          eval_left_shift(dt, e_diff);
 777          res.exponent() = a.exponent() - e_diff;
 778          eval_subtract(dt, b.bits());
 779       }
 780       else if (a.exponent() == (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent() + 1)
 781       {
 782          if (eval_lsb(b.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
 783          {
 784             eval_left_shift(dt, 1);
 785             eval_decrement(dt);
 786             res.exponent() = a.exponent() - 1;
 787          }
 788          else
 789             res.exponent() = a.exponent();
 790       }
 791       else
 792          res.exponent() = a.exponent();
 793    }
 794    else
 795    {
 796       dt = b.bits();
 797       if (b.exponent() <= (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + a.exponent())
 798       {
 799          typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e_diff = a.exponent() - b.exponent();
 800          eval_left_shift(dt, -e_diff);
 801          res.exponent() = b.exponent() + e_diff;
 802          eval_subtract(dt, a.bits());
 803       }
 804       else if (b.exponent() == (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + a.exponent() + 1)
 805       {
 806          if (eval_lsb(a.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
 807          {
 808             eval_left_shift(dt, 1);
 809             eval_decrement(dt);
 810             res.exponent() = b.exponent() - 1;
 811          }
 812          else
 813             res.exponent() = b.exponent();
 814       }
 815       else
 816          res.exponent() = b.exponent();
 817       s = !s;
 818    }
 819
 820    copy_and_round(res, dt);
 821    if (res.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero)
 822       res.sign() = false;
 823    else if (res.sign() != s)
 824       res.negate();
 825    res.check_invariants();
 826 }
 827
 828 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 829 inline void eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b)
 830 {
 831    if (a.sign() == b.sign())
 832       do_eval_add(res, a, b);
 833    else
 834       do_eval_subtract(res, a, b);
 835 }
 836
 837 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 838 inline void eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a)
 839 {
 840    return eval_add(res, res, a);
 841 }
 842
 843 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 844 inline void eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b)
 845 {
 846    if (a.sign() != b.sign())
 847       do_eval_add(res, a, b);
 848    else
 849       do_eval_subtract(res, a, b);
 850 }
 851
 852 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 853 inline void eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a)
 854 {
 855    return eval_subtract(res, res, a);
 856 }
 857
 858 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 859 inline void eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b)
 860 {
 861    using default_ops::eval_bit_test;
 862    using default_ops::eval_multiply;
 863
 864    // Special cases first:
 865    switch (a.exponent())
 866    {
 867    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
 868    {
 869       if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)
 870          res = b;
 871       else if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity)
 872          res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
 873       else
 874       {
 875          bool s     = a.sign() != b.sign();
 876          res        = a;
 877          res.sign() = s;
 878       }
 879       return;
 880    }
 881    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
 882       switch (b.exponent())
 883       {
 884       case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
 885          res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
 886          break;
 887       case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
 888          res = b;
 889          break;
 890       default:
 891          bool s     = a.sign() != b.sign();
 892          res        = a;
 893          res.sign() = s;
 894          break;
 895       }
 896       return;
 897    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
 898       res = a;
 899       return;
 900    }
 901    if (b.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
 902    {
 903       bool s     = a.sign() != b.sign();
 904       res        = b;
 905       res.sign() = s;
 906       return;
 907    }
 908    if ((a.exponent() > 0) && (b.exponent() > 0))
 909    {
 910       if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent + 2 - a.exponent() < b.exponent())
 911       {
 912          // We will certainly overflow:
 913          bool s         = a.sign() != b.sign();
 914          res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
 915          res.sign()     = s;
 916          res.bits()     = static_cast<limb_type>(0u);
 917          return;
 918       }
 919    }
 920    if ((a.exponent() < 0) && (b.exponent() < 0))
 921    {
 922       if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - 2 - a.exponent() > b.exponent())
 923       {
 924          // We will certainly underflow:
 925          res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
 926          res.sign()     = a.sign() != b.sign();
 927          res.bits()     = static_cast<limb_type>(0u);
 928          return;
 929       }
 930    }
 931
 932    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
 933    eval_multiply(dt, a.bits(), b.bits());
 934    res.exponent() = a.exponent() + b.exponent() - (Exponent)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
 935    copy_and_round(res, dt);
 936    res.check_invariants();
 937    res.sign() = a.sign() != b.sign();
 938 }
 939
 940 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 941 inline void eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a)
 942 {
 943    eval_multiply(res, res, a);
 944 }
 945
 946 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
 947 inline typename enable_if_c<is_unsigned<U>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const U& b)
 948 {
 949    using default_ops::eval_bit_test;
 950    using default_ops::eval_multiply;
 951
 952    // Special cases first:
 953    switch (a.exponent())
 954    {
 955    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
 956    {
 957       bool s     = a.sign();
 958       res        = a;
 959       res.sign() = s;
 960       return;
 961    }
 962    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
 963       if (b == 0)
 964          res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
 965       else
 966          res = a;
 967       return;
 968    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
 969       res = a;
 970       return;
 971    }
 972
 973    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type                                                                     dt;
 974    typedef typename boost::multiprecision::detail::canonical<U, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::type canon_ui_type;
 975    eval_multiply(dt, a.bits(), static_cast<canon_ui_type>(b));
 976    res.exponent() = a.exponent();
 977    copy_and_round(res, dt);
 978    res.check_invariants();
 979    res.sign() = a.sign();
 980 }
 981
 982 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
 983 inline typename enable_if_c<is_unsigned<U>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const U& b)
 984 {
 985    eval_multiply(res, res, b);
 986 }
 987
 988 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
 989 inline typename enable_if_c<is_signed<S>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const S& b)
 990 {
 991    typedef typename make_unsigned<S>::type ui_type;
 992    eval_multiply(res, a, static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(b)));
 993    if (b < 0)
 994       res.negate();
 995 }
 996
 997 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
 998 inline typename enable_if_c<is_signed<S>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const S& b)
 999 {
1000    eval_multiply(res, res, b);
1001 }
1002
1003 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1004 inline void eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& u, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& v)
1005 {
1006 #ifdef BOOST_MSVC
1007 #pragma warning(push)
1008 #pragma warning(disable : 6326) // comparison of two constants
1009 #endif
1010    using default_ops::eval_bit_test;
1011    using default_ops::eval_get_sign;
1012    using default_ops::eval_increment;
1013    using default_ops::eval_qr;
1014    using default_ops::eval_subtract;
1015
1016    //
1017    // Special cases first:
1018    //
1019    switch (u.exponent())
1020    {
1021    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1022    {
1023       switch (v.exponent())
1024       {
1025       case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1026       case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1027          res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1028          return;
1029       }
1030       bool s     = u.sign() != v.sign();
1031       res        = u;
1032       res.sign() = s;
1033       return;
1034    }
1035    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1036    {
1037       switch (v.exponent())
1038       {
1039       case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1040       case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1041          res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1042          return;
1043       }
1044       bool s     = u.sign() != v.sign();
1045       res        = u;
1046       res.sign() = s;
1047       return;
1048    }
1049    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1050       res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1051       return;
1052    }
1053    switch (v.exponent())
1054    {
1055    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1056    {
1057       bool s     = u.sign() != v.sign();
1058       res        = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
1059       res.sign() = s;
1060       return;
1061    }
1062    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1063       res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
1064       res.bits()     = limb_type(0);
1065       res.sign()     = u.sign() != v.sign();
1066       return;
1067    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1068       res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1069       return;
1070    }
1071
1072    // We can scale u and v so that both are integers, then perform integer
1073    // division to obtain quotient q and remainder r, such that:
1074    //
1075    // q * v + r = u
1076    //
1077    // and hense:
1078    //
1079    // q + r/v = u/v
1080    //
1081    // From this, assuming q has cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count
1082    // bits we only need to determine whether
1083    // r/v is less than, equal to, or greater than 0.5 to determine rounding -
1084    // this we can do with a shift and comparison.
1085    //
1086    // We can set the exponent and sign of the result up front:
1087    //
1088    if ((v.exponent() < 0) && (u.exponent() > 0))
1089    {
1090       // Check for overflow:
1091       if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent + v.exponent() < u.exponent() - 1)
1092       {
1093          res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
1094          res.sign()     = u.sign() != v.sign();
1095          res.bits()     = static_cast<limb_type>(0u);
1096          return;
1097       }
1098    }
1099    else if ((v.exponent() > 0) && (u.exponent() < 0))
1100    {
1101       // Check for underflow:
1102       if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent + v.exponent() > u.exponent())
1103       {
1104          // We will certainly underflow:
1105          res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
1106          res.sign()     = u.sign() != v.sign();
1107          res.bits()     = static_cast<limb_type>(0u);
1108          return;
1109       }
1110    }
1111    res.exponent() = u.exponent() - v.exponent() - 1;
1112    res.sign()     = u.sign() != v.sign();
1113    //
1114    // Now get the quotient and remainder:
1115    //
1116    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(u.bits()), t2(v.bits()), q, r;
1117    eval_left_shift(t, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count);
1118    eval_qr(t, t2, q, r);
1119    //
1120    // We now have either "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count"
1121    // or "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1" significant
1122    // bits in q.
1123    //
1124    static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
1125    if (eval_bit_test(q, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
1126    {
1127       //
1128       // OK we have cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1 bits,
1129       // so we already have rounding info,
1130       // we just need to changes things if the last bit is 1 and either the
1131       // remainder is non-zero (ie we do not have a tie) or the quotient would
1132       // be odd if it were shifted to the correct number of bits (ie a tiebreak).
1133       //
1134       BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count));
1135       if ((q.limbs()[0] & 1u) && (eval_get_sign(r) || (q.limbs()[0] & 2u)))
1136       {
1137          eval_increment(q);
1138       }
1139    }
1140    else
1141    {
1142       //
1143       // We have exactly "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" bits in q.
1144       // Get rounding info, which we can get by comparing 2r with v.
1145       // We want to call copy_and_round to handle rounding and general cleanup,
1146       // so we'll left shift q and add some fake digits on the end to represent
1147       // how we'll be rounding.
1148       //
1149       BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
1150       static const unsigned lshift = (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count < limb_bits) ? 2 : limb_bits;
1151       eval_left_shift(q, lshift);
1152       res.exponent() -= lshift;
1153       eval_left_shift(r, 1u);
1154       int c = r.compare(v.bits());
1155       if (c == 0)
1156          q.limbs()[0] |= static_cast<limb_type>(1u) << (lshift - 1);
1157       else if (c > 0)
1158          q.limbs()[0] |= (static_cast<limb_type>(1u) << (lshift - 1)) + static_cast<limb_type>(1u);
1159    }
1160    copy_and_round(res, q);
1161 #ifdef BOOST_MSVC
1162 #pragma warning(pop)
1163 #endif
1164 }
1165
1166 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1167 inline void eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1168 {
1169    eval_divide(res, res, arg);
1170 }
1171
1172 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
1173 inline typename enable_if_c<is_unsigned<U>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& u, const U& v)
1174 {
1175 #ifdef BOOST_MSVC
1176 #pragma warning(push)
1177 #pragma warning(disable : 6326) // comparison of two constants
1178 #endif
1179    using default_ops::eval_bit_test;
1180    using default_ops::eval_get_sign;
1181    using default_ops::eval_increment;
1182    using default_ops::eval_qr;
1183    using default_ops::eval_subtract;
1184
1185    //
1186    // Special cases first:
1187    //
1188    switch (u.exponent())
1189    {
1190    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1191    {
1192       if (v == 0)
1193       {
1194          res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1195          return;
1196       }
1197       bool s     = u.sign() != (v < 0);
1198       res        = u;
1199       res.sign() = s;
1200       return;
1201    }
1202    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1203       res = u;
1204       return;
1205    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1206       res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1207       return;
1208    }
1209    if (v == 0)
1210    {
1211       bool s     = u.sign();
1212       res        = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
1213       res.sign() = s;
1214       return;
1215    }
1216
1217    // We can scale u and v so that both are integers, then perform integer
1218    // division to obtain quotient q and remainder r, such that:
1219    //
1220    // q * v + r = u
1221    //
1222    // and hense:
1223    //
1224    // q + r/v = u/v
1225    //
1226    // From this, assuming q has "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count, we only need to determine whether
1227    // r/v is less than, equal to, or greater than 0.5 to determine rounding -
1228    // this we can do with a shift and comparison.
1229    //
1230    // We can set the exponent and sign of the result up front:
1231    //
1232    int gb         = msb(v);
1233    res.exponent() = u.exponent() - static_cast<Exponent>(gb) - static_cast<Exponent>(1);
1234    res.sign()     = u.sign();
1235    //
1236    // Now get the quotient and remainder:
1237    //
1238    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(u.bits()), q, r;
1239    eval_left_shift(t, gb + 1);
1240    eval_qr(t, number<typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::canonical_value(v), q, r);
1241    //
1242    // We now have either "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" or "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1" significant cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in q.
1243    //
1244    static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
1245    if (eval_bit_test(q, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
1246    {
1247       //
1248       // OK we have cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1 cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count, so we already have rounding info,
1249       // we just need to changes things if the last bit is 1 and the
1250       // remainder is non-zero (ie we do not have a tie).
1251       //
1252       BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count));
1253       if ((q.limbs()[0] & 1u) && eval_get_sign(r))
1254       {
1255          eval_increment(q);
1256       }
1257    }
1258    else
1259    {
1260       //
1261       // We have exactly "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in q.
1262       // Get rounding info, which we can get by comparing 2r with v.
1263       // We want to call copy_and_round to handle rounding and general cleanup,
1264       // so we'll left shift q and add some fake cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count on the end to represent
1265       // how we'll be rounding.
1266       //
1267       BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
1268       static const unsigned lshift = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count < limb_bits ? 2 : limb_bits;
1269       eval_left_shift(q, lshift);
1270       res.exponent() -= lshift;
1271       eval_left_shift(r, 1u);
1272       int c = r.compare(number<typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::canonical_value(v));
1273       if (c == 0)
1274          q.limbs()[0] |= static_cast<limb_type>(1u) << (lshift - 1);
1275       else if (c > 0)
1276          q.limbs()[0] |= (static_cast<limb_type>(1u) << (lshift - 1)) + static_cast<limb_type>(1u);
1277    }
1278    copy_and_round(res, q);
1279 #ifdef BOOST_MSVC
1280 #pragma warning(pop)
1281 #endif
1282 }
1283
1284 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
1285 inline typename enable_if_c<is_unsigned<U>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const U& v)
1286 {
1287    eval_divide(res, res, v);
1288 }
1289
1290 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
1291 inline typename enable_if_c<is_signed<S>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& u, const S& v)
1292 {
1293    typedef typename make_unsigned<S>::type ui_type;
1294    eval_divide(res, u, static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(v)));
1295    if (v < 0)
1296       res.negate();
1297 }
1298
1299 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
1300 inline typename enable_if_c<is_signed<S>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const S& v)
1301 {
1302    eval_divide(res, res, v);
1303 }
1304
1305 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1306 inline int eval_get_sign(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1307 {
1308    return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero ? 0 : arg.sign() ? -1 : 1;
1309 }
1310
1311 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1312 inline bool eval_is_zero(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1313 {
1314    return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
1315 }
1316
1317 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1318 inline bool eval_eq(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b)
1319 {
1320    if (a.exponent() == b.exponent())
1321    {
1322       if (a.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero)
1323          return true;
1324       return (a.sign() == b.sign()) && (a.bits().compare(b.bits()) == 0) && (a.exponent() != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan);
1325    }
1326    return false;
1327 }
1328
1329 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1330 inline void eval_convert_to(boost::long_long_type* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1331 {
1332    switch (arg.exponent())
1333    {
1334    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1335       *res = 0;
1336       return;
1337    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1338       BOOST_THROW_EXCEPTION(std::runtime_error("Could not convert NaN to integer."));
1339    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1340       *res = (std::numeric_limits<boost::long_long_type>::max)();
1341       if (arg.sign())
1342          *res = -*res;
1343       return;
1344    }
1345    typedef typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift_type;
1346    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type                                                                                                                                                              man(arg.bits());
1347    shift_type                                                                                                                                                                                                                                        shift = (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - arg.exponent();
1348    if (shift > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
1349    {
1350       *res = 0;
1351       return;
1352    }
1353    if (arg.sign() && (arg.compare((std::numeric_limits<boost::long_long_type>::min)()) <= 0))
1354    {
1355       *res = (std::numeric_limits<boost::long_long_type>::min)();
1356       return;
1357    }
1358    else if (!arg.sign() && (arg.compare((std::numeric_limits<boost::long_long_type>::max)()) >= 0))
1359    {
1360       *res = (std::numeric_limits<boost::long_long_type>::max)();
1361       return;
1362    }
1363
1364    if (shift < 0)
1365    {
1366       if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - shift <= std::numeric_limits<boost::long_long_type>::digits)
1367       {
1368          // We have more bits in long_long_type than the float, so it's OK to left shift:
1369          eval_convert_to(res, man);
1370          *res <<= -shift;
1371       }
1372       else
1373       {
1374          *res = (std::numeric_limits<boost::long_long_type>::max)();
1375          return;
1376       }
1377    }
1378    else
1379    {
1380       eval_right_shift(man, shift);
1381       eval_convert_to(res, man);
1382    }
1383    if (arg.sign())
1384    {
1385       *res = -*res;
1386    }
1387 }
1388
1389 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1390 inline void eval_convert_to(boost::ulong_long_type* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1391 {
1392    switch (arg.exponent())
1393    {
1394    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1395       *res = 0;
1396       return;
1397    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1398       BOOST_THROW_EXCEPTION(std::runtime_error("Could not convert NaN to integer."));
1399    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1400       *res = (std::numeric_limits<boost::ulong_long_type>::max)();
1401       return;
1402    }
1403    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type                                                                                                                                                              man(arg.bits());
1404    typedef typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift_type;
1405    shift_type                                                                                                                                                                                                                                        shift = (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - arg.exponent();
1406    if (shift > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
1407    {
1408       *res = 0;
1409       return;
1410    }
1411    else if (shift < 0)
1412    {
1413       if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - shift <= std::numeric_limits<boost::ulong_long_type>::digits)
1414       {
1415          // We have more bits in ulong_long_type than the float, so it's OK to left shift:
1416          eval_convert_to(res, man);
1417          *res <<= -shift;
1418          return;
1419       }
1420       *res = (std::numeric_limits<boost::ulong_long_type>::max)();
1421       return;
1422    }
1423    eval_right_shift(man, shift);
1424    eval_convert_to(res, man);
1425 }
1426
1427 template <class Float, unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1428 inline typename boost::enable_if_c<boost::is_float<Float>::value>::type eval_convert_to(Float* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& original_arg)
1429 {
1430    typedef cpp_bin_float<std::numeric_limits<Float>::digits, digit_base_2, void, Exponent, MinE, MaxE> conv_type;
1431    typedef typename common_type<typename conv_type::exponent_type, int>::type                          common_exp_type;
1432    //
1433    // Special cases first:
1434    //
1435    switch (original_arg.exponent())
1436    {
1437    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1438       *res = 0;
1439       if (original_arg.sign())
1440          *res = -*res;
1441       return;
1442    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1443       *res = std::numeric_limits<Float>::quiet_NaN();
1444       return;
1445    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1446       *res = (std::numeric_limits<Float>::infinity)();
1447       if (original_arg.sign())
1448          *res = -*res;
1449       return;
1450    }
1451    //
1452    // Check for super large exponent that must be converted to infinity:
1453    //
1454    if (original_arg.exponent() > std::numeric_limits<Float>::max_exponent)
1455    {
1456       *res = std::numeric_limits<Float>::has_infinity ? std::numeric_limits<Float>::infinity() : (std::numeric_limits<Float>::max)();
1457       if (original_arg.sign())
1458          *res = -*res;
1459       return;
1460    }
1461    //
1462    // Figure out how many digits we will have in our result,
1463    // allowing for a possibly denormalized result:
1464    //
1465    common_exp_type digits_to_round_to = std::numeric_limits<Float>::digits;
1466    if (original_arg.exponent() < std::numeric_limits<Float>::min_exponent - 1)
1467    {
1468       common_exp_type diff = original_arg.exponent();
1469       diff -= std::numeric_limits<Float>::min_exponent - 1;
1470       digits_to_round_to += diff;
1471    }
1472    if (digits_to_round_to < 0)
1473    {
1474       // Result must be zero:
1475       *res = 0;
1476       if (original_arg.sign())
1477          *res = -*res;
1478       return;
1479    }
1480    //
1481    // Perform rounding first, then afterwards extract the digits:
1482    //
1483    cpp_bin_float<std::numeric_limits<Float>::digits, digit_base_2, Allocator, Exponent, MinE, MaxE> arg;
1484    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type             bits(original_arg.bits());
1485    arg.exponent() = original_arg.exponent();
1486    copy_and_round(arg, bits, (int)digits_to_round_to);
1487    common_exp_type e = arg.exponent();
1488    e -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1;
1489    static const unsigned limbs_needed      = std::numeric_limits<Float>::digits / (sizeof(*arg.bits().limbs()) * CHAR_BIT) + (std::numeric_limits<Float>::digits % (sizeof(*arg.bits().limbs()) * CHAR_BIT) ? 1 : 0);
1490    unsigned              first_limb_needed = arg.bits().size() - limbs_needed;
1491    *res                                    = 0;
1492    e += first_limb_needed * sizeof(*arg.bits().limbs()) * CHAR_BIT;
1493    while (first_limb_needed < arg.bits().size())
1494    {
1495       *res += std::ldexp(static_cast<Float>(arg.bits().limbs()[first_limb_needed]), static_cast<int>(e));
1496       ++first_limb_needed;
1497       e += sizeof(*arg.bits().limbs()) * CHAR_BIT;
1498    }
1499    if (original_arg.sign())
1500       *res = -*res;
1501 }
1502
1503 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1504 inline void eval_frexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, Exponent* e)
1505 {
1506    switch (arg.exponent())
1507    {
1508    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1509    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1510    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1511       *e  = 0;
1512       res = arg;
1513       return;
1514    }
1515    res            = arg;
1516    *e             = arg.exponent() + 1;
1517    res.exponent() = -1;
1518 }
1519
1520 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
1521 inline void eval_frexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, I* pe)
1522 {
1523    Exponent e;
1524    eval_frexp(res, arg, &e);
1525    if ((e > (std::numeric_limits<I>::max)()) || (e < (std::numeric_limits<I>::min)()))
1526    {
1527       BOOST_THROW_EXCEPTION(std::runtime_error("Exponent was outside of the range of the argument type to frexp."));
1528    }
1529    *pe = static_cast<I>(e);
1530 }
1531
1532 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1533 inline void eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, Exponent e)
1534 {
1535    switch (arg.exponent())
1536    {
1537    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1538    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1539    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1540       res = arg;
1541       return;
1542    }
1543    if ((e > 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent - e < arg.exponent()))
1544    {
1545       // Overflow:
1546       res        = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
1547       res.sign() = arg.sign();
1548    }
1549    else if ((e < 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - e > arg.exponent()))
1550    {
1551       // Underflow:
1552       res = limb_type(0);
1553    }
1554    else
1555    {
1556       res = arg;
1557       res.exponent() += e;
1558    }
1559 }
1560
1561 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
1562 inline typename enable_if_c<is_unsigned<I>::value>::type eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, I e)
1563 {
1564    typedef typename make_signed<I>::type si_type;
1565    if (e > static_cast<I>((std::numeric_limits<si_type>::max)()))
1566       res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
1567    else
1568       eval_ldexp(res, arg, static_cast<si_type>(e));
1569 }
1570
1571 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
1572 inline typename enable_if_c<is_signed<I>::value>::type eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, I e)
1573 {
1574    if ((e > (std::numeric_limits<Exponent>::max)()) || (e < (std::numeric_limits<Exponent>::min)()))
1575    {
1576       res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
1577       if (e < 0)
1578          res.negate();
1579    }
1580    else
1581       eval_ldexp(res, arg, static_cast<Exponent>(e));
1582 }
1583
1584 /*
1585 * Sign manipulation
1586 */
1587
1588 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1589 inline void eval_abs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1590 {
1591    res        = arg;
1592    res.sign() = false;
1593 }
1594
1595 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1596 inline void eval_fabs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1597 {
1598    res        = arg;
1599    res.sign() = false;
1600 }
1601
1602 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1603 inline int eval_fpclassify(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1604 {
1605    switch (arg.exponent())
1606    {
1607    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1608       return FP_ZERO;
1609    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1610       return FP_INFINITE;
1611    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1612       return FP_NAN;
1613    }
1614    return FP_NORMAL;
1615 }
1616
1617 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1618 inline void eval_sqrt(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1619 {
1620    using default_ops::eval_bit_test;
1621    using default_ops::eval_increment;
1622    using default_ops::eval_integer_sqrt;
1623    switch (arg.exponent())
1624    {
1625    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1626       errno = EDOM;
1627       // fallthrough...
1628    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1629       res = arg;
1630       return;
1631    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1632       if (arg.sign())
1633       {
1634          res   = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1635          errno = EDOM;
1636       }
1637       else
1638          res = arg;
1639       return;
1640    }
1641    if (arg.sign())
1642    {
1643       res   = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1644       errno = EDOM;
1645       return;
1646    }
1647
1648    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(arg.bits()), r, s;
1649    eval_left_shift(t, arg.exponent() & 1 ? cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count : cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1);
1650    eval_integer_sqrt(s, r, t);
1651
1652    if (!eval_bit_test(s, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
1653    {
1654       // We have exactly the right number of cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in the result, round as required:
1655       if (s.compare(r) < 0)
1656       {
1657          eval_increment(s);
1658       }
1659    }
1660    typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type ae = arg.exponent();
1661    res.exponent()                                                                               = ae / 2;
1662    if ((ae & 1) && (ae < 0))
1663       --res.exponent();
1664    copy_and_round(res, s);
1665 }
1666
1667 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1668 inline void eval_floor(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1669 {
1670    using default_ops::eval_increment;
1671    switch (arg.exponent())
1672    {
1673    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1674       errno = EDOM;
1675       // fallthrough...
1676    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1677    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1678       res = arg;
1679       return;
1680    }
1681    typedef typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift_type;
1682    shift_type                                                                                                                                                                                                                                        shift =
1683        (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - arg.exponent() - 1;
1684    if ((arg.exponent() > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) || (shift <= 0))
1685    {
1686       // Either arg is already an integer, or a special value:
1687       res = arg;
1688       return;
1689    }
1690    if (shift >= (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
1691    {
1692       res = static_cast<signed_limb_type>(arg.sign() ? -1 : 0);
1693       return;
1694    }
1695    bool fractional = (shift_type)eval_lsb(arg.bits()) < shift;
1696    res             = arg;
1697    eval_right_shift(res.bits(), shift);
1698    if (fractional && res.sign())
1699    {
1700       eval_increment(res.bits());
1701       if (eval_msb(res.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - shift)
1702       {
1703          // Must have extended result by one bit in the increment:
1704          --shift;
1705          ++res.exponent();
1706       }
1707    }
1708    eval_left_shift(res.bits(), shift);
1709 }
1710
1711 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1712 inline void eval_ceil(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
1713 {
1714    using default_ops::eval_increment;
1715    switch (arg.exponent())
1716    {
1717    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1718       errno = EDOM;
1719       // fallthrough...
1720    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1721    case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1722       res = arg;
1723       return;
1724    }
1725    typedef typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift_type;
1726    shift_type                                                                                                                                                                                                                                        shift = (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - arg.exponent() - 1;
1727    if ((arg.exponent() > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) || (shift <= 0))
1728    {
1729       // Either arg is already an integer, or a special value:
1730       res = arg;
1731       return;
1732    }
1733    if (shift >= (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
1734    {
1735       bool s     = arg.sign(); // takes care of signed zeros
1736       res        = static_cast<signed_limb_type>(arg.sign() ? 0 : 1);
1737       res.sign() = s;
1738       return;
1739    }
1740    bool fractional = (shift_type)eval_lsb(arg.bits()) < shift;
1741    res             = arg;
1742    eval_right_shift(res.bits(), shift);
1743    if (fractional && !res.sign())
1744    {
1745       eval_increment(res.bits());
1746       if (eval_msb(res.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - shift)
1747       {
1748          // Must have extended result by one bit in the increment:
1749          --shift;
1750          ++res.exponent();
1751       }
1752    }
1753    eval_left_shift(res.bits(), shift);
1754 }
1755
1756 template <unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2>
1757 int eval_signbit(const cpp_bin_float<D1, B1, A1, E1, M1, M2>& val)
1758 {
1759    return val.sign();
1760 }
1761
1762 template <unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2>
1763 inline std::size_t hash_value(const cpp_bin_float<D1, B1, A1, E1, M1, M2>& val)
1764 {
1765    std::size_t result = hash_value(val.bits());
1766    boost::hash_combine(result, val.exponent());
1767    boost::hash_combine(result, val.sign());
1768    return result;
1769 }
1770
1771 } // namespace backends
1772
1773 #ifdef BOOST_NO_SFINAE_EXPR
1774
1775 namespace detail {
1776
1777 template <unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2, unsigned D2, backends::digit_base_type B2, class A2, class E2, E2 M3, E2 M4>
1778 struct is_explicitly_convertible<backends::cpp_bin_float<D1, B1, A1, E1, M1, M2>, backends::cpp_bin_float<D2, B2, A2, E2, M3, M4> > : public mpl::true_
1779 {};
1780 template <class FloatT, unsigned D2, backends::digit_base_type B2, class A2, class E2, E2 M3, E2 M4>
1781 struct is_explicitly_convertible<FloatT, backends::cpp_bin_float<D2, B2, A2, E2, M3, M4> > : public boost::is_floating_point<FloatT>
1782 {};
1783
1784 } // namespace detail
1785 #endif
1786
1787 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Exponent, Exponent MinE, Exponent MaxE, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
1788 inline boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>
1789     copysign BOOST_PREVENT_MACRO_SUBSTITUTION(
1790         const boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>& a,
1791         const boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>& b)
1792 {
1793    boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> res(a);
1794    res.backend().sign() = b.backend().sign();
1795    return res;
1796 }
1797
1798 using backends::cpp_bin_float;
1799 using backends::digit_base_10;
1800 using backends::digit_base_2;
1801
1802 template <unsigned Digits, backends::digit_base_type DigitBase, class Exponent, Exponent MinE, Exponent MaxE, class Allocator>
1803 struct number_category<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > : public boost::mpl::int_<boost::multiprecision::number_kind_floating_point>
1804 {};
1805
1806 template <unsigned Digits, backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1807 struct expression_template_default<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >
1808 {
1809    static const expression_template_option value = is_void<Allocator>::value ? et_off : et_on;
1810 };
1811
1812 typedef number<backends::cpp_bin_float<50> >  cpp_bin_float_50;
1813 typedef number<backends::cpp_bin_float<100> > cpp_bin_float_100;
1814
1815 typedef number<backends::cpp_bin_float<24, backends::digit_base_2, void, boost::int16_t, -126, 127>, et_off>        cpp_bin_float_single;
1816 typedef number<backends::cpp_bin_float<53, backends::digit_base_2, void, boost::int16_t, -1022, 1023>, et_off>      cpp_bin_float_double;
1817 typedef number<backends::cpp_bin_float<64, backends::digit_base_2, void, boost::int16_t, -16382, 16383>, et_off>    cpp_bin_float_double_extended;
1818 typedef number<backends::cpp_bin_float<113, backends::digit_base_2, void, boost::int16_t, -16382, 16383>, et_off>   cpp_bin_float_quad;
1819 typedef number<backends::cpp_bin_float<237, backends::digit_base_2, void, boost::int32_t, -262142, 262143>, et_off> cpp_bin_float_oct;
1820
1821 } // namespace multiprecision
1822
1823 namespace math {
1824
1825 using boost::multiprecision::copysign;
1826 using boost::multiprecision::signbit;
1827
1828 } // namespace math
1829
1830 } // namespace boost
1831
1832 #include <boost/multiprecision/cpp_bin_float/io.hpp>
1833 #include <boost/multiprecision/cpp_bin_float/transcendental.hpp>
1834
1835 namespace std {
1836
1837 //
1838 // numeric_limits [partial] specializations for the types declared in this header:
1839 //
1840 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1841 class numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >
1842 {
1843    typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> number_type;
1844
1845  public:
1846    BOOST_STATIC_CONSTEXPR bool is_specialized = true;
1847    static number_type(min)()
1848    {
1849       initializer.do_nothing();
1850       static std::pair<bool, number_type> value;
1851       if (!value.first)
1852       {
1853          value.first = true;
1854          typedef typename boost::mpl::front<typename number_type::backend_type::unsigned_types>::type ui_type;
1855          value.second.backend()            = ui_type(1u);
1856          value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
1857       }
1858       return value.second;
1859    }
1860    static number_type(max)()
1861    {
1862       initializer.do_nothing();
1863       static std::pair<bool, number_type> value;
1864       if (!value.first)
1865       {
1866          value.first = true;
1867          if (boost::is_void<Allocator>::value)
1868             eval_complement(value.second.backend().bits(), value.second.backend().bits());
1869          else
1870          {
1871             // We jump through hoops here using the backend type directly just to keep VC12 happy
1872             // (ie compiler workaround, for very strange compiler bug):
1873             using boost::multiprecision::default_ops::eval_add;
1874             using boost::multiprecision::default_ops::eval_decrement;
1875             using boost::multiprecision::default_ops::eval_left_shift;
1876             typedef typename number_type::backend_type::rep_type                                int_backend_type;
1877             typedef typename boost::mpl::front<typename int_backend_type::unsigned_types>::type ui_type;
1878             int_backend_type                                                                    i;
1879             i = ui_type(1u);
1880             eval_left_shift(i, boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1);
1881             int_backend_type j(i);
1882             eval_decrement(i);
1883             eval_add(j, i);
1884             value.second.backend().bits() = j;
1885          }
1886          value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
1887       }
1888       return value.second;
1889    }
1890    BOOST_STATIC_CONSTEXPR number_type lowest()
1891    {
1892       return -(max)();
1893    }
1894    BOOST_STATIC_CONSTEXPR int digits   = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count;
1895    BOOST_STATIC_CONSTEXPR int digits10 = (digits - 1) * 301 / 1000;
1896    // Is this really correct???
1897    BOOST_STATIC_CONSTEXPR int  max_digits10 = (digits * 301 / 1000) + 3;
1898    BOOST_STATIC_CONSTEXPR bool is_signed    = true;
1899    BOOST_STATIC_CONSTEXPR bool is_integer   = false;
1900    BOOST_STATIC_CONSTEXPR bool is_exact     = false;
1901    BOOST_STATIC_CONSTEXPR int  radix        = 2;
1902    static number_type          epsilon()
1903    {
1904       initializer.do_nothing();
1905       static std::pair<bool, number_type> value;
1906       if (!value.first)
1907       {
1908          // We jump through hoops here just to keep VC12 happy (ie compiler workaround, for very strange compiler bug):
1909          typedef typename boost::mpl::front<typename number_type::backend_type::unsigned_types>::type ui_type;
1910          value.first            = true;
1911          value.second.backend() = ui_type(1u);
1912          value.second           = ldexp(value.second, 1 - (int)digits);
1913       }
1914       return value.second;
1915    }
1916    // What value should this be????
1917    static number_type round_error()
1918    {
1919       // returns 0.5
1920       initializer.do_nothing();
1921       static std::pair<bool, number_type> value;
1922       if (!value.first)
1923       {
1924          value.first = true;
1925          // We jump through hoops here just to keep VC12 happy (ie compiler workaround, for very strange compiler bug):
1926          typedef typename boost::mpl::front<typename number_type::backend_type::unsigned_types>::type ui_type;
1927          value.second.backend() = ui_type(1u);
1928          value.second           = ldexp(value.second, -1);
1929       }
1930       return value.second;
1931    }
1932    BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type min_exponent      = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
1933    BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type min_exponent10    = (min_exponent / 1000) * 301L;
1934    BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type max_exponent      = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
1935    BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type max_exponent10    = (max_exponent / 1000) * 301L;
1936    BOOST_STATIC_CONSTEXPR bool                                                                                                             has_infinity      = true;
1937    BOOST_STATIC_CONSTEXPR bool                                                                                                             has_quiet_NaN     = true;
1938    BOOST_STATIC_CONSTEXPR bool                                                                                                             has_signaling_NaN = false;
1939    BOOST_STATIC_CONSTEXPR float_denorm_style has_denorm                                                                                                      = denorm_absent;
1940    BOOST_STATIC_CONSTEXPR bool               has_denorm_loss                                                                                                 = false;
1941    static number_type                        infinity()
1942    {
1943       initializer.do_nothing();
1944       static std::pair<bool, number_type> value;
1945       if (!value.first)
1946       {
1947          value.first                       = true;
1948          value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
1949       }
1950       return value.second;
1951    }
1952    static number_type quiet_NaN()
1953    {
1954       initializer.do_nothing();
1955       static std::pair<bool, number_type> value;
1956       if (!value.first)
1957       {
1958          value.first                       = true;
1959          value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan;
1960       }
1961       return value.second;
1962    }
1963    BOOST_STATIC_CONSTEXPR number_type signaling_NaN()
1964    {
1965       return number_type(0);
1966    }
1967    BOOST_STATIC_CONSTEXPR number_type denorm_min() { return number_type(0); }
1968    BOOST_STATIC_CONSTEXPR bool        is_iec559         = false;
1969    BOOST_STATIC_CONSTEXPR bool        is_bounded        = true;
1970    BOOST_STATIC_CONSTEXPR bool        is_modulo         = false;
1971    BOOST_STATIC_CONSTEXPR bool        traps             = true;
1972    BOOST_STATIC_CONSTEXPR bool        tinyness_before   = false;
1973    BOOST_STATIC_CONSTEXPR float_round_style round_style = round_to_nearest;
1974
1975  private:
1976    struct data_initializer
1977    {
1978       data_initializer()
1979       {
1980          std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::epsilon();
1981          std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::round_error();
1982          (std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::min)();
1983          (std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::max)();
1984          std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity();
1985          std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN();
1986       }
1987       void do_nothing() const {}
1988    };
1989    static const data_initializer initializer;
1990 };
1991
1992 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1993 const typename numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::data_initializer numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::initializer;
1994
1995 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
1996
1997 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1998 BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::digits;
1999 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2000 BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::digits10;
2001 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2002 BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_digits10;
2003 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2004 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_signed;
2005 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2006 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_integer;
2007 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2008 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_exact;
2009 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2010 BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::radix;
2011 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2012 BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::min_exponent;
2013 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2014 BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::min_exponent10;
2015 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2016 BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_exponent;
2017 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2018 BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_exponent10;
2019 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2020 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_infinity;
2021 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2022 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_quiet_NaN;
2023 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2024 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_signaling_NaN;
2025 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2026 BOOST_CONSTEXPR_OR_CONST float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_denorm;
2027 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2028 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_denorm_loss;
2029 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2030 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_iec559;
2031 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2032 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_bounded;
2033 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2034 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_modulo;
2035 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2036 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::traps;
2037 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2038 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::tinyness_before;
2039 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
2040 BOOST_CONSTEXPR_OR_CONST float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::round_style;
2041
2042 #endif
2043
2044 } // namespace std
2045
2046 #endif