133 changes: 119 additions & 14 deletions include/boost/multiprecision/detail/functions/pow.hpp
View file
Expand Up @@ -208,6 +208,7 @@ void eval_exp(T& result, const T& x)
if(type == (int)FP_NAN)
{
result = x;
errno = EDOM;
return;
}
else if(type == (int)FP_INFINITE)
Expand Down Expand Up @@ -326,6 +327,29 @@ void eval_log(T& result, const T& arg)
typedef typename T::exponent_type exp_type;
typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type;
typedef typename mpl::front<typename T::float_types>::type fp_type;
int s = eval_signbit(arg);
switch(eval_fpclassify(arg))
{
case FP_NAN:
result = arg;
errno = EDOM;
return;
case FP_INFINITE:
if(s) break;
result = arg;
return;
case FP_ZERO:
result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend();
result.negate();
errno = ERANGE;
return;
}
if(s)
{
result = std::numeric_limits<number<T> >::quiet_NaN().backend();
errno = EDOM;
return;
}

exp_type e;
T t;
Expand Down Expand Up @@ -427,19 +451,21 @@ inline void eval_pow(T& result, const T& x, const T& a)
return;
}

if(a.compare(si_type(1)) == 0)
if((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0))
{
result = x;
return;
}
if(a.compare(si_type(0)) == 0)
{
result = si_type(1);
return;
}

int type = eval_fpclassify(x);

switch(type)
{
case FP_INFINITE:
result = x;
return;
case FP_ZERO:
switch(eval_fpclassify(a))
{
Expand All @@ -449,13 +475,41 @@ inline void eval_pow(T& result, const T& x, const T& a)
case FP_NAN:
result = a;
break;
case FP_NORMAL:
{
// Need to check for a an odd integer as a special case:
try
{
boost::intmax_t i;
eval_convert_to(&i, a);
if((a.compare(i) == 0) && (i & 1) && eval_signbit(a))
{
result = std::numeric_limits<number<T> >::infinity().backend();
if(eval_signbit(x))
result.negate();
errno = ERANGE;
return;
}
}
catch(const std::exception&)
{
// fallthrough..
}
}
default:
result = x;
if(eval_signbit(a))
{
result = std::numeric_limits<number<T> >::infinity().backend();
errno = ERANGE;
}
else
result = x;
break;
}
return;
case FP_NAN:
result = x;
errno = ERANGE;
return;
default: ;
}
Expand All @@ -478,6 +532,16 @@ inline void eval_pow(T& result, const T& x, const T& a)
}

typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type an;
typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type max_an =
std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ?
(std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::max)() :
static_cast<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>(1) << (sizeof(typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type) * CHAR_BIT - 2);
typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type min_an =
std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ?
(std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::min)() :
-min_an;


T fa;
#ifndef BOOST_NO_EXCEPTIONS
try
Expand Down Expand Up @@ -521,8 +585,33 @@ inline void eval_pow(T& result, const T& x, const T& a)
// conversion failed, just fall through, value is not an integer.
}
#endif
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
eval_floor(result, a);
// -1^INF is a special case in C99:
if((x.compare(si_type(-1)) == 0) && (eval_fpclassify(a) == FP_INFINITE))
{
result = si_type(1);
}
else if(a.compare(result) == 0)
{
// exponent is so large we have no fractional part:
if(x.compare(si_type(-1)) < 0)
{
result = std::numeric_limits<number<T, et_on> >::infinity().backend();
}
else
{
result = si_type(0);
}
}
else if(type == FP_INFINITE)
{
result = std::numeric_limits<number<T, et_on> >::infinity().backend();
}
else if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
{
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
errno = EDOM;
}
else
{
BOOST_THROW_EXCEPTION(std::domain_error("Result of pow is undefined or non-real and there is no NaN for this number type."));
Expand All @@ -534,7 +623,7 @@ inline void eval_pow(T& result, const T& x, const T& a)

eval_subtract(da, a, an);

if((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0))
if((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an))
{
if(a.compare(fp_type(1e-5f)) <= 0)
{
Expand Down Expand Up @@ -618,14 +707,20 @@ void eval_exp2(T& result, const T& arg)
// Check for pure-integer arguments which can be either signed or unsigned.
typename boost::multiprecision::detail::canonical<typename T::exponent_type, T>::type i;
T temp;
eval_trunc(temp, arg);
eval_convert_to(&i, temp);
if(arg.compare(i) == 0)
{
temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
eval_ldexp(result, temp, i);
return;
try {
eval_trunc(temp, arg);
eval_convert_to(&i, temp);
if(arg.compare(i) == 0)
{
temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
eval_ldexp(result, temp, i);
return;
}
}
catch(const boost::math::rounding_error&)
{ /* Fallthrough */ }
catch(const std::runtime_error&)
{ /* Fallthrough */ }

temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(2u);
eval_pow(result, temp, arg);
Expand Down Expand Up @@ -669,6 +764,8 @@ namespace detail{
switch(eval_fpclassify(x))
{
case FP_NAN:
errno = EDOM;
// fallthrough...
case FP_INFINITE:
if(p_sinh)
*p_sinh = x;
Expand Down Expand Up @@ -742,6 +839,14 @@ inline void eval_tanh(T& result, const T& x)
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tanh function is only valid for floating point types.");
T c;
detail::sinhcosh(x, &result, &c);
if((eval_fpclassify(result) == FP_INFINITE) && (eval_fpclassify(c) == FP_INFINITE))
{
bool s = eval_signbit(result) != eval_signbit(c);
result = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
if(s)
result.negate();
return;
}
eval_divide(result, c);
}

Expand Down
57 changes: 47 additions & 10 deletions include/boost/multiprecision/detail/functions/trig.hpp
View file
Expand Up @@ -91,12 +91,15 @@ void eval_sin(T& result, const T& x)
case FP_INFINITE:
case FP_NAN:
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
{
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
errno = EDOM;
}
else
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
return;
case FP_ZERO:
result = ui_type(0);
result = x;
return;
default: ;
}
Expand Down Expand Up @@ -238,7 +241,10 @@ void eval_cos(T& result, const T& x)
case FP_INFINITE:
case FP_NAN:
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
{
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
errno = EDOM;
}
else
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
return;
Expand Down Expand Up @@ -422,12 +428,15 @@ void eval_asin(T& result, const T& x)
case FP_NAN:
case FP_INFINITE:
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
{
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
errno = EDOM;
}
else
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
return;
case FP_ZERO:
result = ui_type(0);
result = x;
return;
default: ;
}
Expand All @@ -442,7 +451,10 @@ void eval_asin(T& result, const T& x)
if(c > 0)
{
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
{
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
errno = EDOM;
}
else
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
return;
Expand Down Expand Up @@ -535,7 +547,10 @@ inline void eval_acos(T& result, const T& x)
case FP_NAN:
case FP_INFINITE:
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
{
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
errno = EDOM;
}
else
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
return;
Expand All @@ -551,7 +566,10 @@ inline void eval_acos(T& result, const T& x)
if(c > 0)
{
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
{
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
errno = EDOM;
}
else
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
return;
Expand Down Expand Up @@ -584,9 +602,10 @@ void eval_atan(T& result, const T& x)
{
case FP_NAN:
result = x;
errno = EDOM;
return;
case FP_ZERO:
result = ui_type(0);
result = x;
return;
case FP_INFINITE:
if(eval_get_sign(x) < 0)
Expand Down Expand Up @@ -694,24 +713,41 @@ void eval_atan2(T& result, const T& y, const T& x)
{
case FP_NAN:
result = y;
errno = EDOM;
return;
case FP_ZERO:
{
int c = eval_get_sign(x);
if(c < 0)
if(eval_signbit(x))
{
result = get_constant_pi<T>();
else if(c >= 0)
result = ui_type(0); // Note we allow atan2(0,0) to be zero, even though it's mathematically undefined
if(eval_signbit(y))
result.negate();
}
else
{
result = y; // Note we allow atan2(0,0) to be +-zero, even though it's mathematically undefined
}
return;
}
case FP_INFINITE:
{
if(eval_fpclassify(x) == FP_INFINITE)
{
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
if(eval_signbit(x))
{
// 3Pi/4
eval_ldexp(result, get_constant_pi<T>(), -2);
eval_subtract(result, get_constant_pi<T>());
if(eval_get_sign(y) >= 0)
result.negate();
}
else
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
{
// Pi/4
eval_ldexp(result, get_constant_pi<T>(), -2);
if(eval_get_sign(y) < 0)
result.negate();
}
}
else
{
Expand All @@ -727,6 +763,7 @@ void eval_atan2(T& result, const T& y, const T& x)
{
case FP_NAN:
result = x;
errno = EDOM;
return;
case FP_ZERO:
{
Expand Down
20 changes: 5 additions & 15 deletions include/boost/multiprecision/float128.hpp
View file
Expand Up @@ -399,15 +399,6 @@ inline void eval_fabs(float128_backend& result, const float128_backend& arg)

inline void eval_trunc(float128_backend& result, const float128_backend& arg)
{
if(isnanq(arg.value()) || isinfq(arg.value()))
{
result = boost::math::policies::raise_rounding_error(
"boost::multiprecision::trunc<%1%>(%1%)", 0,
number<float128_backend, et_off>(arg),
number<float128_backend, et_off>(arg),
boost::math::policies::policy<>()).backend();
return;
}
result.value() = truncq(arg.value());
}
/*
Expand Down Expand Up @@ -494,6 +485,11 @@ inline void eval_multiply_add(float128_backend& result, const float128_backend&
result.value() = fmaq(a.value(), b.value(), c.value());
}

inline int eval_signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const float128_backend& arg)
{
return ::signbitq(arg.value());
}

inline std::size_t hash_value(const float128_backend& val)
{
return boost::hash_value(static_cast<double>(val.value()));
Expand Down Expand Up @@ -552,12 +548,6 @@ inline std::size_t hash_value(const float128_backend& val)
return log1pq(arg.backend().value());
}

template <multiprecision::expression_template_option ExpressionTemplates>
inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates>& arg)
{
return ::signbitq(arg.backend().value());
}

template <multiprecision::expression_template_option ExpressionTemplates>
inline boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates>& a, const boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates>& b)
{
Expand Down
7 changes: 7 additions & 0 deletions include/boost/multiprecision/logged_adaptor.hpp
View file
Expand Up @@ -506,6 +506,13 @@ NON_MEMBER_OP3(fmod, "fmod");
NON_MEMBER_OP3(pow, "pow");
NON_MEMBER_OP3(atan2, "atan2");

template <class Backend>
int eval_signbit(const logged_adaptor<Backend>& val)
{
using default_ops::eval_signbit;
return eval_signbit(val.value());
}

template <class Backend>
std::size_t hash_value(const logged_adaptor<Backend>& val)
{
Expand Down
23 changes: 6 additions & 17 deletions include/boost/multiprecision/mpfr.hpp
View file
Expand Up @@ -1286,11 +1286,6 @@ inline void eval_floor(mpfr_float_backend<Digits10, AllocateType>& result, const
template <unsigned Digits10, mpfr_allocation_type AllocateType>
inline void eval_trunc(mpfr_float_backend<Digits10, AllocateType>& result, const mpfr_float_backend<Digits10, AllocateType>& val)
{
if(0 == mpfr_number_p(val.data()))
{
result = boost::math::policies::raise_rounding_error("boost::multiprecision::trunc<%1%>(%1%)", 0, number<mpfr_float_backend<Digits10, AllocateType> >(val), number<mpfr_float_backend<Digits10, AllocateType> >(val), boost::math::policies::policy<>()).backend();
return;
}
mpfr_trunc(result.data(), val.data());
}
template <unsigned Digits10, mpfr_allocation_type AllocateType>
Expand Down Expand Up @@ -1504,6 +1499,12 @@ inline void eval_multiply_subtract(mpfr_float_backend<Digits10, AllocateType>& r
mpfr_fms(result.data(), a.data(), b.data(), c.data(), GMP_RNDN);
}

template <unsigned Digits10, mpfr_allocation_type AllocateType>
inline int eval_signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const mpfr_float_backend<Digits10, AllocateType>& arg)
{
return (arg.data()[0]._mpfr_sign < 0) ? 1 : 0;
}

template <unsigned Digits10, mpfr_allocation_type AllocateType>
inline std::size_t hash_value(const mpfr_float_backend<Digits10, AllocateType>& val)
{
Expand Down Expand Up @@ -1545,24 +1546,12 @@ typedef number<mpfr_float_backend<0> > mpfr_float;
typedef number<mpfr_float_backend<50, allocate_stack> > static_mpfr_float_50;
typedef number<mpfr_float_backend<100, allocate_stack> > static_mpfr_float_100;

template<unsigned Digits10, boost::multiprecision::mpfr_allocation_type AllocateType, boost::multiprecision::expression_template_option ExpressionTemplates>
inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates>& arg)
{
return (arg.backend().data()[0]._mpfr_sign < 0) ? 1 : 0;
}

template<unsigned Digits10, boost::multiprecision::mpfr_allocation_type AllocateType, boost::multiprecision::expression_template_option ExpressionTemplates>
inline boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates>& a, const boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates>& b)
{
return (boost::multiprecision::signbit)(a) != (boost::multiprecision::signbit)(b) ? boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType>, ExpressionTemplates>(-a) : a;
}

template<unsigned Digits10, boost::multiprecision::mpfr_allocation_type AllocateType, boost::multiprecision::expression_template_option ExpressionTemplates>
inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType> >, ExpressionTemplates>& arg)
{
return (arg.backend().value().data()[0]._mpfr_sign < 0) ? 1 : 0;
}

template<unsigned Digits10, boost::multiprecision::mpfr_allocation_type AllocateType, boost::multiprecision::expression_template_option ExpressionTemplates>
inline boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType> >, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType> >, ExpressionTemplates>& a, const boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpfr_float_backend<Digits10, AllocateType> >, ExpressionTemplates>& b)
{
Expand Down
9 changes: 5 additions & 4 deletions test/test_round.cpp
View file
Expand Up @@ -340,7 +340,7 @@ void test()
#endif
if(std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_THROW(static_cast<T>(round(std::numeric_limits<T>::infinity())), boost::math::rounding_error);
BOOST_CHECK_EQUAL(static_cast<T>(round(std::numeric_limits<T>::infinity())), std::numeric_limits<T>::infinity()); // See C99 Annex F.
BOOST_CHECK_THROW(static_cast<T>(iround(std::numeric_limits<T>::infinity())), boost::math::rounding_error);
BOOST_CHECK_THROW(static_cast<T>(iround(-std::numeric_limits<T>::infinity())), boost::math::rounding_error);
BOOST_CHECK_THROW(static_cast<T>(lround(std::numeric_limits<T>::infinity())), boost::math::rounding_error);
Expand All @@ -352,7 +352,7 @@ void test()
}
if(std::numeric_limits<T>::has_quiet_NaN)
{
BOOST_CHECK_THROW(static_cast<T>(round(std::numeric_limits<T>::quiet_NaN())), boost::math::rounding_error);
BOOST_CHECK((boost::multiprecision::isnan)(round(std::numeric_limits<T>::quiet_NaN())));
BOOST_CHECK_THROW(static_cast<T>(iround(std::numeric_limits<T>::quiet_NaN())), boost::math::rounding_error);
BOOST_CHECK_THROW(static_cast<T>(lround(std::numeric_limits<T>::quiet_NaN())), boost::math::rounding_error);
#ifdef BOOST_HAS_LONG_LONG
Expand All @@ -369,7 +369,8 @@ void test()
#endif
if(std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_THROW(static_cast<T>(trunc(std::numeric_limits<T>::infinity())), boost::math::rounding_error);
BOOST_CHECK_EQUAL(static_cast<T>(trunc(std::numeric_limits<T>::infinity())), std::numeric_limits<T>::infinity());
BOOST_CHECK_EQUAL(static_cast<T>(trunc(-std::numeric_limits<T>::infinity())), -std::numeric_limits<T>::infinity());
BOOST_CHECK_THROW(static_cast<T>(itrunc(std::numeric_limits<T>::infinity())), boost::math::rounding_error);
BOOST_CHECK_THROW(static_cast<T>(itrunc(-std::numeric_limits<T>::infinity())), boost::math::rounding_error);
BOOST_CHECK_THROW(static_cast<T>(ltrunc(std::numeric_limits<T>::infinity())), boost::math::rounding_error);
Expand All @@ -381,7 +382,7 @@ void test()
}
if(std::numeric_limits<T>::has_quiet_NaN)
{
BOOST_CHECK_THROW(static_cast<T>(trunc(std::numeric_limits<T>::quiet_NaN())), boost::math::rounding_error);
BOOST_CHECK((boost::multiprecision::isnan)(trunc(std::numeric_limits<T>::quiet_NaN())));
BOOST_CHECK_THROW(static_cast<T>(itrunc(std::numeric_limits<T>::quiet_NaN())), boost::math::rounding_error);
BOOST_CHECK_THROW(static_cast<T>(ltrunc(std::numeric_limits<T>::quiet_NaN())), boost::math::rounding_error);
#ifdef BOOST_HAS_LONG_LONG
Expand Down
1,410 changes: 1,410 additions & 0 deletions test/test_sf_import_c99.cpp
View file

Large diffs are not rendered by default.