/usr/include/boost/multiprecision
Edit: /usr/include/boost/multiprecision/cpp_bin_float.hpp (98102B)
///////////////////////////////////////////////////////////////
// Copyright 2013 John Maddock. Distributed under the Boost
// Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
#ifndef BOOST_MATH_CPP_BIN_FLOAT_HPP
#define BOOST_MATH_CPP_BIN_FLOAT_HPP
#include
#include
#include
#include
//
// Some includes we need from Boost.Math, since we rely on that library to provide these functions:
//
#include
#include
#include
#include
#include
#include
#ifdef BOOST_HAS_FLOAT128
#include
#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
inline typename enable_if_c::value, bool>::type is_negative(U) { return false; }
template
inline typename disable_if_c::value, bool>::type is_negative(S s) { return s < 0; }
template ::value == number_kind_floating_point>
struct is_cpp_bin_float_implicitly_constructible_from_type
{
static const bool value = false;
};
template
struct is_cpp_bin_float_implicitly_constructible_from_type
{
static const bool value = (std::numeric_limits::digits <= (int)bit_count) && (std::numeric_limits::radix == 2) && std::numeric_limits::is_specialized
#ifdef BOOST_HAS_FLOAT128
&& !boost::is_same::value
#endif
&& (is_floating_point::value || is_number::value);
};
template ::value == number_kind_floating_point>
struct is_cpp_bin_float_explicitly_constructible_from_type
{
static const bool value = false;
};
template
struct is_cpp_bin_float_explicitly_constructible_from_type
{
static const bool value = (std::numeric_limits::digits > (int)bit_count) && (std::numeric_limits::radix == 2) && std::numeric_limits::is_specialized
#ifdef BOOST_HAS_FLOAT128
&& !boost::is_same::value
#endif
;
};
} // namespace detail
template
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::value ? bit_count : 0, bit_count, is_void::value ? unsigned_magnitude : signed_magnitude, unchecked, Allocator> rep_type;
typedef cpp_int_backend::value ? 2 * bit_count : 0, 2 * bit_count, is_void::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_types;
typedef Exponent exponent_type;
static const exponent_type max_exponent_limit = boost::integer_traits::const_max - 2 * static_cast(bit_count);
static const exponent_type min_exponent_limit = boost::integer_traits::const_min + 2 * static_cast(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())))
: m_data(o.m_data), m_exponent(o.m_exponent), m_sign(o.m_sign) {}
template
cpp_bin_float(const cpp_bin_float& o, typename boost::enable_if_c<(bit_count >= cpp_bin_float::bit_count)>::type const* = 0)
{
*this = o;
}
template
explicit cpp_bin_float(const cpp_bin_float& o, typename boost::disable_if_c<(bit_count >= cpp_bin_float::bit_count)>::type const* = 0)
: m_exponent(o.exponent()), m_sign(o.sign())
{
*this = o;
}
template
cpp_bin_float(const Float& f,
typename boost::enable_if_c::value>::type const* = 0)
: m_data(), m_exponent(0), m_sign(false)
{
this->assign_float(f);
}
template
explicit cpp_bin_float(const Float& f,
typename boost::enable_if_c::value>::type const* = 0)
: m_data(), m_exponent(0), m_sign(false)
{
this->assign_float(f);
}
#ifdef BOOST_HAS_FLOAT128
template
cpp_bin_float(const Float& f,
typename boost::enable_if_c<
boost::is_same::value && ((int)bit_count >= 113)>::type const* = 0)
: m_data(), m_exponent(0), m_sign(false)
{
this->assign_float(f);
}
template
explicit cpp_bin_float(const Float& f,
typename boost::enable_if_c<
boost::is_same::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() = std::declval()))
{
m_data = o.m_data;
m_exponent = o.m_exponent;
m_sign = o.m_sign;
return *this;
}
template
cpp_bin_float& operator=(const cpp_bin_float& 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::rep_type b(f.bits());
this->exponent() = f.exponent() + (E)bit_count - (E)cpp_bin_float::bit_count;
this->sign() = f.sign();
copy_and_round(*this, b);
}
return *this;
}
#ifdef BOOST_HAS_FLOAT128
template
typename boost::enable_if_c<
(number_category::value == number_kind_floating_point)
//&& (std::numeric_limits::digits <= (int)bit_count)
&& ((std::numeric_limits::radix == 2) || (boost::is_same::value)),
cpp_bin_float&>::type
operator=(const Float& f)
#else
template
typename boost::enable_if_c<
(number_category::value == number_kind_floating_point)
//&& (std::numeric_limits::digits <= (int)bit_count)
&& (std::numeric_limits::radix == 2),
cpp_bin_float&>::type
operator=(const Float& f)
#endif
{
return assign_float(f);
}
#ifdef BOOST_HAS_FLOAT128
template
typename boost::enable_if_c::value, cpp_bin_float&>::type assign_float(Float f)
{
using default_ops::eval_add;
typedef typename boost::multiprecision::detail::canonical::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::type ui_type;
m_data = static_cast(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(ipart);
eval_add(*this, t);
}
m_exponent += static_cast(e);
return *this;
}
#endif
#ifdef BOOST_HAS_FLOAT128
template
typename boost::enable_if_c::value && !is_same::value, cpp_bin_float&>::type assign_float(Float f)
#else
template
typename boost::enable_if_c::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::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::type ui_type;
m_data = static_cast(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(f);
#endif
f -= ipart;
m_exponent += bits;
cpp_bin_float t;
t = static_cast(ipart);
eval_add(*this, t);
}
m_exponent += static_cast(e);
return *this;
}
template
typename boost::enable_if_c<
(number_category::value == number_kind_floating_point) && !boost::is_floating_point::value && (number_category::value == number_kind_floating_point),
cpp_bin_float&>::type
assign_float(Float f)
{
BOOST_MATH_STD_USING
using default_ops::eval_add;
using default_ops::eval_convert_to;
using default_ops::eval_get_sign;
using default_ops::eval_subtract;
typedef typename boost::multiprecision::detail::canonical::type f_int_type;
typedef typename boost::multiprecision::detail::canonical::type bf_int_type;
switch (eval_fpclassify(f))
{
case FP_ZERO:
m_data = limb_type(0);
m_sign = (eval_get_sign(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 = eval_get_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::type ui_type;
m_data = static_cast(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(ipart));
m_exponent += bits;
eval_add(*this, static_cast(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 = (eval_get_sign(f) > 0);
}
else if (eval_get_sign(m_data) == 0)
{
m_exponent = exponent_zero;
m_sign = (eval_get_sign(f) > 0);
}
return *this;
}
template
cpp_bin_float& assign_float(const number& f)
{
return assign_float(f.backend());
}
template
typename boost::enable_if, cpp_bin_float&>::type operator=(const I& i)
{
using default_ops::eval_bit_test;
if (!i)
{
m_data = static_cast(0);
m_exponent = exponent_zero;
m_sign = false;
}
else
{
typedef typename make_unsigned::type ui_type;
ui_type fi = static_cast(boost::multiprecision::detail::unsigned_abs(i));
typedef typename boost::multiprecision::detail::canonical::type ar_type;
m_data = static_cast(fi);
unsigned shift = msb(fi);
if (shift >= bit_count)
{
m_exponent = static_cast(shift);
m_data = static_cast(fi >> (shift + 1 - bit_count));
}
else
{
m_exponent = static_cast(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
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
void serialize(Archive& ar, const unsigned int /*version*/)
{
ar& boost::make_nvp("data", m_data);
ar& boost::make_nvp("exponent", m_exponent);
ar& boost::make_nvp("sign", m_sign);
}
};
#ifdef BOOST_MSVC
#pragma warning(pop)
#endif
template
inline void copy_and_round(cpp_bin_float& res, Int& arg, int bits_to_keep = cpp_bin_float::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_bit_test;
using default_ops::eval_get_sign;
using default_ops::eval_increment;
using default_ops::eval_left_shift;
using default_ops::eval_lsb;
using default_ops::eval_msb;
using default_ops::eval_right_shift;
// cancellation may have resulted in arg being all zeros:
if (eval_get_sign(arg) == 0)
{
res.exponent() = cpp_bin_float::exponent_zero;
res.sign() = false;
res.bits() = static_cast(0u);
return;
}
int msb = eval_msb(arg);
if (static_cast(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(bits_to_keep - msb - 1);
}
else if (static_cast(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(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::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::bit_count - bits_to_keep);
res.exponent() -= static_cast(cpp_bin_float::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::exponent_zero;
return;
}
// Result must be normalized:
BOOST_ASSERT(((int)eval_msb(res.bits()) == cpp_bin_float::bit_count - 1));
if (res.exponent() > cpp_bin_float::max_exponent)
{
// Overflow:
res.exponent() = cpp_bin_float::exponent_infinity;
res.bits() = static_cast(0u);
}
else if (res.exponent() < cpp_bin_float::min_exponent)
{
// Underflow:
res.exponent() = cpp_bin_float::exponent_zero;
res.bits() = static_cast(0u);
}
}
template
inline void do_eval_add(cpp_bin_float& res, const cpp_bin_float& a, const cpp_bin_float& 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::exponent_type exponent_type;
typename cpp_bin_float::double_rep_type dt;
// Special cases first:
switch (a.exponent())
{
case cpp_bin_float::exponent_zero:
{
bool s = a.sign();
res = b;
res.sign() = s;
return;
}
case cpp_bin_float::exponent_infinity:
if (b.exponent() == cpp_bin_float::exponent_nan)
res = b;
else
res = a;
return; // result is still infinite.
case cpp_bin_float::exponent_nan:
res = a;
return; // result is still a NaN.
}
switch (b.exponent())
{
case cpp_bin_float::exponent_zero:
res = a;
return;
case cpp_bin_float::exponent_infinity:
res = b;
if (res.sign())
res.negate();
return; // result is infinite.
case cpp_bin_float::exponent_nan:
res = b;
return; // result is a NaN.
}
BOOST_STATIC_ASSERT(boost::integer_traits::const_max - cpp_bin_float::bit_count > cpp_bin_float::max_exponent);
bool s = a.sign();
dt = a.bits();
if (a.exponent() > (int)cpp_bin_float::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
inline void do_eval_subtract(cpp_bin_float& res, const cpp_bin_float& a, const cpp_bin_float& b)
{
using default_ops::eval_bit_test;
using default_ops::eval_decrement;
using default_ops::eval_subtract;
typename cpp_bin_float::double_rep_type dt;
// Special cases first:
switch (a.exponent())
{
case cpp_bin_float::exponent_zero:
if (b.exponent() == cpp_bin_float::exponent_nan)
res = std::numeric_limits > >::quiet_NaN().backend();
else
{
bool s = a.sign();
res = b;
if (res.exponent() == cpp_bin_float::exponent_zero)
res.sign() = false;
else if (res.sign() == s)
res.negate();
}
return;
case cpp_bin_float::exponent_infinity:
if ((b.exponent() == cpp_bin_float::exponent_nan) || (b.exponent() == cpp_bin_float::exponent_infinity))
res = std::numeric_limits > >::quiet_NaN().backend();
else
res = a;
return;
case cpp_bin_float::exponent_nan:
res = a;
return; // result is still a NaN.
}
switch (b.exponent())
{
case cpp_bin_float::exponent_zero:
res = a;
return;
case cpp_bin_float::exponent_infinity:
res.exponent() = cpp_bin_float::exponent_infinity;
res.sign() = !a.sign();
res.bits() = static_cast(0u);
return; // result is a NaN.
case cpp_bin_float::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::bit_count + b.exponent())
{
typename cpp_bin_float::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::bit_count + b.exponent() + 1)
{
if (eval_lsb(b.bits()) != cpp_bin_float::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::bit_count + a.exponent())
{
typename cpp_bin_float::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::bit_count + a.exponent() + 1)
{
if (eval_lsb(a.bits()) != cpp_bin_float::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::exponent_zero)
res.sign() = false;
else if (res.sign() != s)
res.negate();
res.check_invariants();
}
template
inline void eval_add(cpp_bin_float& res, const cpp_bin_float& a, const cpp_bin_float& b)
{
if (a.sign() == b.sign())
do_eval_add(res, a, b);
else
do_eval_subtract(res, a, b);
}
template
inline void eval_add(cpp_bin_float& res, const cpp_bin_float& a)
{
return eval_add(res, res, a);
}
template
inline void eval_subtract(cpp_bin_float& res, const cpp_bin_float& a, const cpp_bin_float& b)
{
if (a.sign() != b.sign())
do_eval_add(res, a, b);
else
do_eval_subtract(res, a, b);
}
template
inline void eval_subtract(cpp_bin_float& res, const cpp_bin_float& a)
{
return eval_subtract(res, res, a);
}
template
inline void eval_multiply(cpp_bin_float& res, const cpp_bin_float& a, const cpp_bin_float& b)
{
using default_ops::eval_bit_test;
using default_ops::eval_multiply;
// Special cases first:
switch (a.exponent())
{
case cpp_bin_float::exponent_zero:
{
if (b.exponent() == cpp_bin_float::exponent_nan)
res = b;
else if (b.exponent() == cpp_bin_float::exponent_infinity)
res = std::numeric_limits > >::quiet_NaN().backend();
else
{
bool s = a.sign() != b.sign();
res = a;
res.sign() = s;
}
return;
}
case cpp_bin_float::exponent_infinity:
switch (b.exponent())
{
case cpp_bin_float::exponent_zero:
res = std::numeric_limits > >::quiet_NaN().backend();
break;
case cpp_bin_float::exponent_nan:
res = b;
break;
default:
bool s = a.sign() != b.sign();
res = a;
res.sign() = s;
break;
}
return;
case cpp_bin_float::exponent_nan:
res = a;
return;
}
if (b.exponent() > cpp_bin_float::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::max_exponent + 2 - a.exponent() < b.exponent())
{
// We will certainly overflow:
bool s = a.sign() != b.sign();
res.exponent() = cpp_bin_float::exponent_infinity;
res.sign() = s;
res.bits() = static_cast(0u);
return;
}
}
if ((a.exponent() < 0) && (b.exponent() < 0))
{
if (cpp_bin_float::min_exponent - 2 - a.exponent() > b.exponent())
{
// We will certainly underflow:
res.exponent() = cpp_bin_float::exponent_zero;
res.sign() = a.sign() != b.sign();
res.bits() = static_cast(0u);
return;
}
}
typename cpp_bin_float::double_rep_type dt;
eval_multiply(dt, a.bits(), b.bits());
res.exponent() = a.exponent() + b.exponent() - (Exponent)cpp_bin_float::bit_count + 1;
copy_and_round(res, dt);
res.check_invariants();
res.sign() = a.sign() != b.sign();
}
template
inline void eval_multiply(cpp_bin_float& res, const cpp_bin_float& a)
{
eval_multiply(res, res, a);
}
template
inline typename enable_if_c::value>::type eval_multiply(cpp_bin_float& res, const cpp_bin_float& 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::exponent_zero:
{
bool s = a.sign();
res = a;
res.sign() = s;
return;
}
case cpp_bin_float::exponent_infinity:
if (b == 0)
res = std::numeric_limits > >::quiet_NaN().backend();
else
res = a;
return;
case cpp_bin_float::exponent_nan:
res = a;
return;
}
typename cpp_bin_float::double_rep_type dt;
typedef typename boost::multiprecision::detail::canonical::double_rep_type>::type canon_ui_type;
eval_multiply(dt, a.bits(), static_cast(b));
res.exponent() = a.exponent();
copy_and_round(res, dt);
res.check_invariants();
res.sign() = a.sign();
}
template
inline typename enable_if_c::value>::type eval_multiply(cpp_bin_float& res, const U& b)
{
eval_multiply(res, res, b);
}
template
inline typename enable_if_c::value>::type eval_multiply(cpp_bin_float& res, const cpp_bin_float& a, const S& b)
{
typedef typename make_unsigned::type ui_type;
eval_multiply(res, a, static_cast(boost::multiprecision::detail::unsigned_abs(b)));
if (b < 0)
res.negate();
}
template
inline typename enable_if_c::value>::type eval_multiply(cpp_bin_float& res, const S& b)
{
eval_multiply(res, res, b);
}
template
inline void eval_divide(cpp_bin_float& res, const cpp_bin_float& u, const cpp_bin_float& v)
{
#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable : 6326) // comparison of two constants
#endif
using default_ops::eval_bit_test;
using default_ops::eval_get_sign;
using default_ops::eval_increment;
using default_ops::eval_qr;
using default_ops::eval_subtract;
//
// Special cases first:
//
switch (u.exponent())
{
case cpp_bin_float::exponent_zero:
{
switch (v.exponent())
{
case cpp_bin_float::exponent_zero:
case cpp_bin_float::exponent_nan:
res = std::numeric_limits > >::quiet_NaN().backend();
return;
}
bool s = u.sign() != v.sign();
res = u;
res.sign() = s;
return;
}
case cpp_bin_float::exponent_infinity:
{
switch (v.exponent())
{
case cpp_bin_float::exponent_infinity:
case cpp_bin_float::exponent_nan:
res = std::numeric_limits