///////////////////////////////////////////////////////////////

// Copyright 2013 John Maddock. Distributed under the Boost

// Software License, Version 1.0. (See accompanying file

// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_

#ifndef BOOST_MATH_CPP_BIN_FLOAT_HPP

#define BOOST_MATH_CPP_BIN_FLOAT_HPP

#include <boost/multiprecision/cpp_int.hpp>

#include <boost/multiprecision/integer.hpp>

#include <boost/math/special_functions/trunc.hpp>

#include <boost/multiprecision/detail/float_string_cvt.hpp>

//

// Some includes we need from Boost.Math, since we rely on that library to provide these functions:

//

#include <boost/math/special_functions/asinh.hpp>

#include <boost/math/special_functions/acosh.hpp>

#include <boost/math/special_functions/atanh.hpp>

#include <boost/math/special_functions/cbrt.hpp>

#include <boost/math/special_functions/expm1.hpp>

#include <boost/math/special_functions/gamma.hpp>

#ifdef BOOST_HAS_FLOAT128

#include <quadmath.h>

#endif

namespace boost{ namespace multiprecision{ namespace backends{

enum digit_base_type

{

digit_base_2 = 2,

digit_base_10 = 10

};

#ifdef BOOST_MSVC

#pragma warning(push)

#pragma warning(disable:4522 6326) // multiple assignment operators specified, comparison of two constants

#endif

namespace detail{

template <class U>

inline typename enable_if_c<is_unsigned<U>::value, bool>::type is_negative(U) { return false; }

template <class S>

inline typename disable_if_c<is_unsigned<S>::value, bool>::type is_negative(S s) { return s < 0; }

template <class Float, int, bool = number_category<Float>::value == number_kind_floating_point>

struct is_cpp_bin_float_implicitly_constructible_from_type

{

static const bool value = false;

};

template <class Float, int bit_count>

struct is_cpp_bin_float_implicitly_constructible_from_type<Float, bit_count, true>

{

static const bool value = (std::numeric_limits<Float>::digits <= (int)bit_count)

&& (std::numeric_limits<Float>::radix == 2)

&& std::numeric_limits<Float>::is_specialized

#ifdef BOOST_HAS_FLOAT128

&& !boost::is_same<Float, __float128>::value

#endif

&& (is_floating_point<Float>::value || is_number<Float>::value)

;

};

template <class Float, int, bool = number_category<Float>::value == number_kind_floating_point>

struct is_cpp_bin_float_explicitly_constructible_from_type

{

static const bool value = false;

};

template <class Float, int bit_count>

struct is_cpp_bin_float_explicitly_constructible_from_type<Float, bit_count, true>

{

static const bool value = (std::numeric_limits<Float>::digits > (int)bit_count)

&& (std::numeric_limits<Float>::radix == 2)

&& std::numeric_limits<Float>::is_specialized

#ifdef BOOST_HAS_FLOAT128

&& !boost::is_same<Float, __float128>::value

#endif

;

};

}

template <unsigned Digits, digit_base_type DigitBase = digit_base_10, class Allocator = void, class Exponent = int, Exponent MinExponent = 0, Exponent MaxExponent = 0>

class cpp_bin_float

{

public:

static const unsigned bit_count = DigitBase == digit_base_2 ? Digits : (Digits * 1000uL) / 301uL + (((Digits * 1000uL) % 301) ? 2u : 1u);

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;

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;

typedef typename rep_type::signed_types signed_types;

typedef typename rep_type::unsigned_types unsigned_types;

typedef boost::mpl::list<float, double, long double> float_types;

typedef Exponent exponent_type;

static const exponent_type max_exponent_limit = boost::integer_traits<exponent_type>::const_max - 2 * static_cast<exponent_type>(bit_count);

static const exponent_type min_exponent_limit = boost::integer_traits<exponent_type>::const_min + 2 * static_cast<exponent_type>(bit_count);

BOOST_STATIC_ASSERT_MSG(MinExponent >= min_exponent_limit, "Template parameter MinExponent is too negative for our internal logic to function correctly, sorry!");

BOOST_STATIC_ASSERT_MSG(MaxExponent <= max_exponent_limit, "Template parameter MaxExponent is too large for our internal logic to function correctly, sorry!");

BOOST_STATIC_ASSERT_MSG(MinExponent <= 0, "Template parameter MinExponent can not be positive!");

BOOST_STATIC_ASSERT_MSG(MaxExponent >= 0, "Template parameter MaxExponent can not be negative!");

static const exponent_type max_exponent = MaxExponent == 0 ? max_exponent_limit : MaxExponent;

static const exponent_type min_exponent = MinExponent == 0 ? min_exponent_limit : MinExponent;

static const exponent_type exponent_zero = max_exponent + 1;

static const exponent_type exponent_infinity = max_exponent + 2;

static const exponent_type exponent_nan = max_exponent + 3;

private:

rep_type m_data;

exponent_type m_exponent;

bool m_sign;

public:

cpp_bin_float() BOOST_MP_NOEXCEPT_IF(noexcept(rep_type())) : m_data(), m_exponent(exponent_zero), m_sign(false) {}

cpp_bin_float(const cpp_bin_float &o) BOOST_MP_NOEXCEPT_IF(noexcept(rep_type(std::declval<const rep_type&>())))

: m_data(o.m_data), m_exponent(o.m_exponent), m_sign(o.m_sign) {}

template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>

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)

{

*this = o;

}

template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>

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)

: m_exponent(o.exponent()), m_sign(o.sign())

{

*this = o;

}

template <class Float>

cpp_bin_float(const Float& f,

typename boost::enable_if_c<detail::is_cpp_bin_float_implicitly_constructible_from_type<Float, bit_count>::value>::type const* = 0)

: m_data(), m_exponent(0), m_sign(false)

{

this->assign_float(f);

}

template <class Float>

explicit cpp_bin_float(const Float& f,

typename boost::enable_if_c<detail::is_cpp_bin_float_explicitly_constructible_from_type<Float, bit_count>::value>::type const* = 0)

: m_data(), m_exponent(0), m_sign(false)

{

this->assign_float(f);

}

#ifdef BOOST_HAS_FLOAT128

template <class Float>

cpp_bin_float(const Float& f,

typename boost::enable_if_c<

boost::is_same<Float, __float128>::value

&& ((int)bit_count >= 113)

>::type const* = 0)

: m_data(), m_exponent(0), m_sign(false)

{

this->assign_float(f);

}

template <class Float>

explicit cpp_bin_float(const Float& f,

typename boost::enable_if_c<

boost::is_same<Float, __float128>::value

&& ((int)bit_count < 113)

>::type const* = 0)

: m_data(), m_exponent(0), m_sign(false)

{

this->assign_float(f);

}

#endif

cpp_bin_float& operator=(const cpp_bin_float &o) BOOST_MP_NOEXCEPT_IF(noexcept(std::declval<rep_type&>() = std::declval<const rep_type&>()))

{

m_data = o.m_data;

m_exponent = o.m_exponent;

m_sign = o.m_sign;

return *this;

}

template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>

cpp_bin_float& operator=(const cpp_bin_float<D, B, A, E, MinE, MaxE> &f)

{

switch(eval_fpclassify(f))

{

case FP_ZERO:

m_data = limb_type(0);

m_sign = f.sign();

m_exponent = exponent_zero;

break;

case FP_NAN:

m_data = limb_type(0);

m_sign = false;

m_exponent = exponent_nan;

break;;

case FP_INFINITE:

m_data = limb_type(0);

m_sign = f.sign();

m_exponent = exponent_infinity;

break;

default:

typename cpp_bin_float<D, B, A, E, MinE, MaxE>::rep_type b(f.bits());

this->exponent() = f.exponent() + (E)bit_count - (E)cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count;

this->sign() = f.sign();

copy_and_round(*this, b);

}

return *this;

}

#ifdef BOOST_HAS_FLOAT128

template <class Float>

typename boost::enable_if_c<

(number_category<Float>::value == number_kind_floating_point)

//&& (std::numeric_limits<Float>::digits <= (int)bit_count)

&& ((std::numeric_limits<Float>::radix == 2) || (boost::is_same<Float, __float128>::value)), cpp_bin_float&>::type

operator=(const Float& f)

#else

template <class Float>

typename boost::enable_if_c<

(number_category<Float>::value == number_kind_floating_point)

//&& (std::numeric_limits<Float>::digits <= (int)bit_count)

&& (std::numeric_limits<Float>::radix == 2), cpp_bin_float&>::type

operator=(const Float& f)

#endif

{

return assign_float(f);

}

#ifdef BOOST_HAS_FLOAT128

template <class Float>

typename boost::enable_if_c<boost::is_same<Float, __float128>::value, cpp_bin_float& >::type assign_float(Float f)

{

using default_ops::eval_add;

typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type;

if(f == 0)

{

m_data = limb_type(0);

m_sign = (signbitq(f) > 0);

m_exponent = exponent_zero;

return *this;

}

else if(isnanq(f))

{

m_data = limb_type(0);

m_sign = false;

m_exponent = exponent_nan;

return *this;

}

else if(isinfq(f))

{

m_data = limb_type(0);

m_sign = (f < 0);

m_exponent = exponent_infinity;

return *this;

}

if(f < 0)

{

*this = -f;

this->negate();

return *this;

}

typedef typename mpl::front<unsigned_types>::type ui_type;

m_data = static_cast<ui_type>(0u);

m_sign = false;

m_exponent = 0;

static const int bits = sizeof(int) * CHAR_BIT - 1;

int e;

f = frexpq(f, &e);

while(f)

{

f = ldexpq(f, bits);

e -= bits;

int ipart = (int)truncq(f);

f -= ipart;

m_exponent += bits;

cpp_bin_float t;

t = static_cast<bf_int_type>(ipart);

eval_add(*this, t);

}

m_exponent += static_cast<Exponent>(e);

return *this;

}

#endif

#ifdef BOOST_HAS_FLOAT128

template <class Float>

typename boost::enable_if_c<is_floating_point<Float>::value && !is_same<Float, __float128>::value, cpp_bin_float&>::type assign_float(Float f)

#else

template <class Float>

typename boost::enable_if_c<is_floating_point<Float>::value, cpp_bin_float&>::type assign_float(Float f)

#endif

{

BOOST_MATH_STD_USING

using default_ops::eval_add;

typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type;

switch((boost::math::fpclassify)(f))

{

case FP_ZERO:

m_data = limb_type(0);

m_sign = ((boost::math::signbit)(f) > 0);

m_exponent = exponent_zero;

return *this;

case FP_NAN:

m_data = limb_type(0);

m_sign = false;

m_exponent = exponent_nan;

return *this;

case FP_INFINITE:

m_data = limb_type(0);

m_sign = (f < 0);

m_exponent = exponent_infinity;

return *this;

}

if(f < 0)

{

*this = -f;

this->negate();

return *this;

}

typedef typename mpl::front<unsigned_types>::type ui_type;

m_data = static_cast<ui_type>(0u);

m_sign = false;

m_exponent = 0;

static const int bits = sizeof(int) * CHAR_BIT - 1;

int e;

f = frexp(f, &e);

while(f)

{

f = ldexp(f, bits);

e -= bits;

#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS

int ipart = itrunc(f);

#else

int ipart = static_cast<int>(f);

#endif

f -= ipart;

m_exponent += bits;

cpp_bin_float t;

t = static_cast<bf_int_type>(ipart);

eval_add(*this, t);

}

m_exponent += static_cast<Exponent>(e);

return *this;

}

template <class Float>

typename boost::enable_if_c<

(number_category<Float>::value == number_kind_floating_point)

&& !boost::is_floating_point<Float>::value

&& is_number<Float>::value,

cpp_bin_float&>::type assign_float(Float f)

{

BOOST_MATH_STD_USING

using default_ops::eval_add;

using default_ops::eval_get_sign;

using default_ops::eval_convert_to;

using default_ops::eval_subtract;

typedef typename boost::multiprecision::detail::canonical<int, Float>::type f_int_type;

typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type;

switch(eval_fpclassify(f))

{

case FP_ZERO:

m_data = limb_type(0);

m_sign = ((boost::math::signbit)(f) > 0);

m_exponent = exponent_zero;

return *this;

case FP_NAN:

m_data = limb_type(0);

m_sign = false;

m_exponent = exponent_nan;

return *this;

case FP_INFINITE:

m_data = limb_type(0);

m_sign = (f < 0);

m_exponent = exponent_infinity;

return *this;

}

if(eval_get_sign(f) < 0)

{

f.negate();

*this = f;

this->negate();

return *this;

}

typedef typename mpl::front<unsigned_types>::type ui_type;

m_data = static_cast<ui_type>(0u);

m_sign = false;

m_exponent = 0;

static const int bits = sizeof(int) * CHAR_BIT - 1;

int e;

eval_frexp(f, f, &e);

while(eval_get_sign(f) != 0)

{

eval_ldexp(f, f, bits);

e -= bits;

int ipart;

eval_convert_to(&ipart, f);

eval_subtract(f, static_cast<f_int_type>(ipart));

m_exponent += bits;

eval_add(*this, static_cast<bf_int_type>(ipart));

}

m_exponent += e;

if(m_exponent > max_exponent)

m_exponent = exponent_infinity;

if(m_exponent < min_exponent)

{

m_data = limb_type(0u);

m_exponent = exponent_zero;

m_sign = ((boost::math::signbit)(f) > 0);

}

else if(eval_get_sign(m_data) == 0)

{

m_exponent = exponent_zero;

m_sign = ((boost::math::signbit)(f) > 0);

}

return *this;

}

template <class I>

typename boost::enable_if<is_integral<I>, cpp_bin_float&>::type operator=(const I& i)

{

using default_ops::eval_bit_test;

if(!i)

{

m_data = static_cast<limb_type>(0);

m_exponent = exponent_zero;

m_sign = false;

}

else

{

typedef typename make_unsigned<I>::type ui_type;

ui_type fi = static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(i));

typedef typename boost::multiprecision::detail::canonical<ui_type, rep_type>::type ar_type;

m_data = static_cast<ar_type>(fi);

unsigned shift = msb(fi);

if(shift >= bit_count)

{

m_exponent = static_cast<Exponent>(shift);

m_data = static_cast<ar_type>(fi >> (shift + 1 - bit_count));

}

else

{

m_exponent = static_cast<Exponent>(shift);

eval_left_shift(m_data, bit_count - shift - 1);

}

BOOST_ASSERT(eval_bit_test(m_data, bit_count-1));

m_sign = detail::is_negative(i);

}

return *this;

}

cpp_bin_float& operator=(const char *s);

void swap(cpp_bin_float &o) BOOST_NOEXCEPT

{

m_data.swap(o.m_data);

std::swap(m_exponent, o.m_exponent);

std::swap(m_sign, o.m_sign);

}

std::string str(std::streamsize dig, std::ios_base::fmtflags f) const;

void negate()

{

if(m_exponent != exponent_nan)

m_sign = !m_sign;

}

int compare(const cpp_bin_float &o) const BOOST_NOEXCEPT

{

if(m_sign != o.m_sign)

return (m_exponent == exponent_zero) && (m_exponent == o.m_exponent) ? 0 : m_sign ? -1 : 1;

int result;

if(m_exponent == exponent_nan)

return -1;

else if(m_exponent != o.m_exponent)

{

if(m_exponent == exponent_zero)

result = -1;

else if(o.m_exponent == exponent_zero)

result = 1;

else

result = m_exponent > o.m_exponent ? 1 : -1;

}

else

result = m_data.compare(o.m_data);

if(m_sign)

result = -result;

return result;

}

template <class A>

int compare(const A& o) const BOOST_NOEXCEPT

{

cpp_bin_float b;

b = o;

return compare(b);

}

rep_type& bits() { return m_data; }

const rep_type& bits()const { return m_data; }

exponent_type& exponent() { return m_exponent; }

const exponent_type& exponent()const { return m_exponent; }

bool& sign() { return m_sign; }

const bool& sign()const { return m_sign; }

void check_invariants()

{

using default_ops::eval_bit_test;

using default_ops::eval_is_zero;

if((m_exponent <= max_exponent) && (m_exponent >= min_exponent))

{

BOOST_ASSERT(eval_bit_test(m_data, bit_count - 1));

}

else

{

BOOST_ASSERT(m_exponent > max_exponent);

BOOST_ASSERT(m_exponent <= exponent_nan);

BOOST_ASSERT(eval_is_zero(m_data));

}

}

template<class Archive>

void serialize(Archive & ar, const unsigned int /*version*/)

{

ar & boost::serialization::make_nvp("data", m_data);

ar & boost::serialization::make_nvp("exponent", m_exponent);

ar & boost::serialization::make_nvp("sign", m_sign);

}

};

#ifdef BOOST_MSVC

#pragma warning(pop)

#endif

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class Int>

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)

{

// Precondition: exponent of res must have been set before this function is called

// as we may need to adjust it based on how many bits_to_keep in arg are set.

using default_ops::eval_msb;

using default_ops::eval_lsb;

using default_ops::eval_left_shift;

using default_ops::eval_bit_test;

using default_ops::eval_right_shift;

using default_ops::eval_increment;

using default_ops::eval_get_sign;

// cancellation may have resulted in arg being all zeros:

if(eval_get_sign(arg) == 0)

{

res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;

res.sign() = false;

res.bits() = static_cast<limb_type>(0u);

return;

}

int msb = eval_msb(arg);

if(static_cast<int>(bits_to_keep) > msb + 1)

{

// Must have had cancellation in subtraction,

// or be converting from a narrower type, so shift left:

res.bits() = arg;

eval_left_shift(res.bits(), bits_to_keep - msb - 1);

res.exponent() -= static_cast<Exponent>(bits_to_keep - msb - 1);

}

else if(static_cast<int>(bits_to_keep) < msb + 1)

{

// We have more bits_to_keep than we need, so round as required,

// first get the rounding bit:

bool roundup = eval_bit_test(arg, msb - bits_to_keep);

// Then check for a tie:

if(roundup && (msb - bits_to_keep == (int)eval_lsb(arg)))

{

// Ties round towards even:

if(!eval_bit_test(arg, msb - bits_to_keep + 1))

roundup = false;

}

// Shift off the bits_to_keep we don't need:

eval_right_shift(arg, msb - bits_to_keep + 1);

res.exponent() += static_cast<Exponent>(msb - bits_to_keep + 1);

if(roundup)

{

eval_increment(arg);

if(bits_to_keep)

{

if(eval_bit_test(arg, bits_to_keep))

{

// This happens very very rairly, all the bits left after

// truncation must be 1's and we're rounding up an order of magnitude:

eval_right_shift(arg, 1u);

++res.exponent();

}

}

else

{

// We get here when bits_to_keep is zero but we're rounding up,

// as a result we end up with a single digit that is a 1:

++bits_to_keep;

}

}

if(bits_to_keep != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)

{

// Normalize result when we're rounding to fewer bits than we can hold, only happens in conversions

// to narrower types:

eval_left_shift(arg, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - bits_to_keep);

res.exponent() -= static_cast<Exponent>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - bits_to_keep);

}

res.bits() = arg;

}

else

{

res.bits() = arg;

}

if(!bits_to_keep && !res.bits().limbs()[0])

{

// We're keeping zero bits and did not round up, so result is zero:

res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;

return;

}

// Result must be normalized:

BOOST_ASSERT(((int)eval_msb(res.bits()) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));

if(res.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)

{

// Overflow:

res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;

res.bits() = static_cast<limb_type>(0u);

}

else if(res.exponent() < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)

{

// Underflow:

res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;

res.bits() = static_cast<limb_type>(0u);

}

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>

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)

{

if(a.exponent() < b.exponent())

{

bool s = a.sign();

do_eval_add(res, b, a);

if(res.sign() != s)

res.negate();

return;

}

using default_ops::eval_add;

using default_ops::eval_bit_test;

typedef typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type exponent_type;

typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;

// Special cases first:

switch(a.exponent())

{

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:

{

bool s = a.sign();

res = b;

res.sign() = s;

return;

}

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:

if(b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)

res = b;

else

res = a;

return; // result is still infinite.

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:

res = a;

return; // result is still a NaN.

}

switch(b.exponent())

{

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:

res = a;

return;

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:

res = b;

if(res.sign())

res.negate();

return; // result is infinite.

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:

res = b;

return; // result is a NaN.

}

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);

bool s = a.sign();

dt = a.bits();

if(a.exponent() > (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent())

{

res.exponent() = a.exponent();

}

else

{

exponent_type e_diff = a.exponent() - b.exponent();

BOOST_ASSERT(e_diff >= 0);

eval_left_shift(dt, e_diff);

res.exponent() = a.exponent() - e_diff;

eval_add(dt, b.bits());

}

copy_and_round(res, dt);

res.check_invariants();

if(res.sign() != s)

res.negate();

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>

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)

{

using default_ops::eval_subtract;

using default_ops::eval_bit_test;

using default_ops::eval_decrement;

typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;

// Special cases first:

switch(a.exponent())

{

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:

if(b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)

res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();

else

{

bool s = a.sign();

res = b;

if(res.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero)

res.sign() = false;

else if(res.sign() == s)

res.negate();

}

return;

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:

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))

res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();

else

res = a;

return;

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:

res = a;

return; // result is still a NaN.

}

switch(b.exponent())

{

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:

res = a;

return;

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:

res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;

res.sign() = !a.sign();

res.bits() = static_cast<limb_type>(0u);

return; // result is a NaN.

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:

res = b;

return; // result is still a NaN.

}

bool s = a.sign();

if((a.exponent() > b.exponent()) || ((a.exponent() == b.exponent()) && a.bits().compare(b.bits()) >= 0))

{

dt = a.bits();

if(a.exponent() <= (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent())

{

typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e_diff = a.exponent() - b.exponent();

eval_left_shift(dt, e_diff);

res.exponent() = a.exponent() - e_diff;

eval_subtract(dt, b.bits());

}

else if(a.exponent() == (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent() + 1)

{

if(eval_lsb(b.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)

{

eval_left_shift(dt, 1);

eval_decrement(dt);

res.exponent() = a.exponent() - 1;

}

else

res.exponent() = a.exponent();

}

else

res.exponent() = a.exponent();

}

else

{

dt = b.bits();

if(b.exponent() <= (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + a.exponent())

{

typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e_diff = a.exponent() - b.exponent();

eval_left_shift(dt, -e_diff);

res.exponent() = b.exponent() + e_diff;

eval_subtract(dt, a.bits());

}

else if(b.exponent() == (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + a.exponent() + 1)

{

if(eval_lsb(a.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)

{

eval_left_shift(dt, 1);

eval_decrement(dt);

res.exponent() = b.exponent() - 1;

}

else

res.exponent() = b.exponent();

}

else

res.exponent() = b.exponent();

s = !s;

}

copy_and_round(res, dt);

if(res.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero)

res.sign() = false;

else if(res.sign() != s)

res.negate();

res.check_invariants();

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>

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)

{

if(a.sign() == b.sign())

do_eval_add(res, a, b);

else

do_eval_subtract(res, a, b);

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>

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)

{

return eval_add(res, res, a);

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>

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)

{

if(a.sign() != b.sign())

do_eval_add(res, a, b);

else

do_eval_subtract(res, a, b);

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>

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)

{

return eval_subtract(res, res, a);

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>

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)

{

using default_ops::eval_bit_test;

using default_ops::eval_multiply;

// Special cases first:

switch(a.exponent())

{

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:

{

if(b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)

res = b;

else if(b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity)

res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();

else

{

bool s = a.sign() != b.sign();

res = a;

res.sign() = s;

}

return;

}

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:

switch(b.exponent())

{

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:

res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();

break;

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:

res = b;

break;

default:

bool s = a.sign() != b.sign();

res = a;

res.sign() = s;

break;

}

return;

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:

res = a;

return;

}

if(b.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)

{

bool s = a.sign() != b.sign();

res = b;

res.sign() = s;

return;

}

if((a.exponent() > 0) && (b.exponent() > 0))

{

if(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent + 2 - a.exponent() < b.exponent())

{

// We will certainly overflow:

bool s = a.sign() != b.sign();

res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;

res.sign() = s;

res.bits() = static_cast<limb_type>(0u);

return;

}

}

if((a.exponent() < 0) && (b.exponent() < 0))

{

if(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - 2 - a.exponent() > b.exponent())

{

// We will certainly underflow:

res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;

res.sign() = a.sign() != b.sign();

res.bits() = static_cast<limb_type>(0u);

return;

}

}

typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;

eval_multiply(dt, a.bits(), b.bits());

res.exponent() = a.exponent() + b.exponent() - (Exponent)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;

copy_and_round(res, dt);

res.check_invariants();

res.sign() = a.sign() != b.sign();

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>

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)

{

eval_multiply(res, res, a);

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>

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)

{

using default_ops::eval_bit_test;

using default_ops::eval_multiply;

// Special cases first:

switch(a.exponent())

{

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:

{

bool s = a.sign();

res = a;

res.sign() = s;

return;

}

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:

if(b == 0)

res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();

else

res = a;

return;

case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:

res = a;

return;

}

typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;

typedef typename boost::multiprecision::detail::canonical<U, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::type canon_ui_type;

eval_multiply(dt, a.bits(), static_cast<canon_ui_type>(b));

res.exponent() = a.exponent();

copy_and_round(res, dt);

res.check_invariants();

res.sign() = a.sign();

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>

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)

{

eval_multiply(res, res, b);

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>

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)

{

typedef typename make_unsigned<S>::type ui_type;

eval_multiply(res, a, static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(b)));

if(b < 0)

res.negate();

}

template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>

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)

{

eval_multiply(res, res, b);