Skip to content

fix for wrong rounding of boost::math::round in classes of corner cases #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Aug 25, 2014
Merged

fix for wrong rounding of boost::math::round in classes of corner cases #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2b5ef4a
pull out for loop body of rounding test for later test
robin-zimmeck Aug 23, 2014
81cfe69
add explicit test cases for rounding corner cases
robin-zimmeck Aug 23, 2014
fb835eb
fix rounding errors
robin-zimmeck Aug 23, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 21 additions & 1 deletion include/boost/math/special_functions/round.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,27 @@ inline typename tools::promote_args<T>::type round(const T& v, const Policy& pol
typedef typename tools::promote_args<T>::type result_type;
if(!(boost::math::isfinite)(v))
return policies::raise_rounding_error("boost::math::round<%1%>(%1%)", 0, static_cast<result_type>(v), static_cast<result_type>(v), pol);
return v < 0 ? static_cast<result_type>(ceil(v - 0.5f)) : static_cast<result_type>(floor(v + 0.5f));

if (-0.5 < v && v < 0.5)
{
// special case to avoid rounding error on the direct
// predecessor of +0.5 resp. the direct successor of -0.5 in
// IEEE floating point types
return 0;
}
else if (v > 0)
{
// subtract v from ceil(v) first in order to avoid rounding
// errors on largest representable integer numbers
result_type c(ceil(v));
return 0.5 < c - v ? c - 1 : c;
}
else
{
// see former branch
result_type f(floor(v));
return 0.5 < v - f ? f + 1 : f;
}
}
template <class T, class Policy>
inline typename tools::promote_args<T>::type round(const T& v, const Policy&, const mpl::true_)
Expand Down
211 changes: 139 additions & 72 deletions test/test_round.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,86 @@ void check_modf_result(T a, T fract, U ipart)
}
}

template <class T>
void test_round_number(T arg)
{
BOOST_MATH_STD_USING
#ifdef BOOST_HAS_LONG_LONG
using boost::math::llround; using boost::math::lltrunc;
#endif

T r = round(arg);
check_within_half(arg, r);
r = trunc(arg);
check_trunc_result(arg, r);
T frac = modf(arg, &r);
check_modf_result(arg, frac, r);

if(abs(r) < (std::numeric_limits<int>::max)())
{
int i = iround(arg);
check_within_half(arg, i);
i = itrunc(arg);
check_trunc_result(arg, i);
r = modf(arg, &i);
check_modf_result(arg, r, i);
}
if(std::numeric_limits<T>::digits >= std::numeric_limits<int>::digits)
{
int si = iround(static_cast<T>((std::numeric_limits<int>::max)()));
check_within_half(static_cast<T>((std::numeric_limits<int>::max)()), si);
si = iround(static_cast<T>((std::numeric_limits<int>::min)()));
check_within_half(static_cast<T>((std::numeric_limits<int>::min)()), si);
si = itrunc(static_cast<T>((std::numeric_limits<int>::max)()));
check_trunc_result(static_cast<T>((std::numeric_limits<int>::max)()), si);
si = itrunc(static_cast<T>((std::numeric_limits<int>::min)()));
check_trunc_result(static_cast<T>((std::numeric_limits<int>::min)()), si);
}
if(abs(r) < (std::numeric_limits<long>::max)())
{
long l = lround(arg);
check_within_half(arg, l);
l = ltrunc(arg);
check_trunc_result(arg, l);
r = modf(arg, &l);
check_modf_result(arg, r, l);
}
if(std::numeric_limits<T>::digits >= std::numeric_limits<long>::digits)
{
long k = lround(static_cast<T>((std::numeric_limits<long>::max)()));
check_within_half(static_cast<T>((std::numeric_limits<long>::max)()), k);
k = lround(static_cast<T>((std::numeric_limits<long>::min)()));
check_within_half(static_cast<T>((std::numeric_limits<long>::min)()), k);
k = ltrunc(static_cast<T>((std::numeric_limits<long>::max)()));
check_trunc_result(static_cast<T>((std::numeric_limits<long>::max)()), k);
k = ltrunc(static_cast<T>((std::numeric_limits<long>::min)()));
check_trunc_result(static_cast<T>((std::numeric_limits<long>::min)()), k);
}

#ifdef BOOST_HAS_LONG_LONG
if(abs(r) < (std::numeric_limits<boost::long_long_type>::max)())
{
boost::long_long_type ll = llround(arg);
check_within_half(arg, ll);
ll = lltrunc(arg);
check_trunc_result(arg, ll);
r = modf(arg, &ll);
check_modf_result(arg, r, ll);
}
if(std::numeric_limits<T>::digits >= std::numeric_limits<boost::long_long_type>::digits)
{
boost::long_long_type j = llround(static_cast<T>((std::numeric_limits<boost::long_long_type>::max)()));
check_within_half(static_cast<T>((std::numeric_limits<boost::long_long_type>::max)()), j);
j = llround(static_cast<T>((std::numeric_limits<boost::long_long_type>::min)()));
check_within_half(static_cast<T>((std::numeric_limits<boost::long_long_type>::min)()), j);
j = lltrunc(static_cast<T>((std::numeric_limits<boost::long_long_type>::max)()));
check_trunc_result(static_cast<T>((std::numeric_limits<boost::long_long_type>::max)()), j);
j = lltrunc(static_cast<T>((std::numeric_limits<boost::long_long_type>::min)()));
check_trunc_result(static_cast<T>((std::numeric_limits<boost::long_long_type>::min)()), j);
}
#endif
}

template <class T>
void test_round(T, const char* name )
{
Expand All @@ -135,76 +215,7 @@ void test_round(T, const char* name )
for(int i = 0; i < 1000; ++i)
{
T arg = get_random<T>();
T r = round(arg);
check_within_half(arg, r);
r = trunc(arg);
check_trunc_result(arg, r);
T frac = modf(arg, &r);
check_modf_result(arg, frac, r);

if(abs(r) < (std::numeric_limits<int>::max)())
{
int i = iround(arg);
check_within_half(arg, i);
i = itrunc(arg);
check_trunc_result(arg, i);
r = modf(arg, &i);
check_modf_result(arg, r, i);
}
if(std::numeric_limits<T>::digits >= std::numeric_limits<int>::digits)
{
int si = iround(static_cast<T>((std::numeric_limits<int>::max)()));
check_within_half(static_cast<T>((std::numeric_limits<int>::max)()), si);
si = iround(static_cast<T>((std::numeric_limits<int>::min)()));
check_within_half(static_cast<T>((std::numeric_limits<int>::min)()), si);
si = itrunc(static_cast<T>((std::numeric_limits<int>::max)()));
check_trunc_result(static_cast<T>((std::numeric_limits<int>::max)()), si);
si = itrunc(static_cast<T>((std::numeric_limits<int>::min)()));
check_trunc_result(static_cast<T>((std::numeric_limits<int>::min)()), si);
}
if(abs(r) < (std::numeric_limits<long>::max)())
{
long l = lround(arg);
check_within_half(arg, l);
l = ltrunc(arg);
check_trunc_result(arg, l);
r = modf(arg, &l);
check_modf_result(arg, r, l);
}
if(std::numeric_limits<T>::digits >= std::numeric_limits<long>::digits)
{
long k = lround(static_cast<T>((std::numeric_limits<long>::max)()));
check_within_half(static_cast<T>((std::numeric_limits<long>::max)()), k);
k = lround(static_cast<T>((std::numeric_limits<long>::min)()));
check_within_half(static_cast<T>((std::numeric_limits<long>::min)()), k);
k = ltrunc(static_cast<T>((std::numeric_limits<long>::max)()));
check_trunc_result(static_cast<T>((std::numeric_limits<long>::max)()), k);
k = ltrunc(static_cast<T>((std::numeric_limits<long>::min)()));
check_trunc_result(static_cast<T>((std::numeric_limits<long>::min)()), k);
}

#ifdef BOOST_HAS_LONG_LONG
if(abs(r) < (std::numeric_limits<boost::long_long_type>::max)())
{
boost::long_long_type ll = llround(arg);
check_within_half(arg, ll);
ll = lltrunc(arg);
check_trunc_result(arg, ll);
r = modf(arg, &ll);
check_modf_result(arg, r, ll);
}
if(std::numeric_limits<T>::digits >= std::numeric_limits<boost::long_long_type>::digits)
{
boost::long_long_type j = llround(static_cast<T>((std::numeric_limits<boost::long_long_type>::max)()));
check_within_half(static_cast<T>((std::numeric_limits<boost::long_long_type>::max)()), j);
j = llround(static_cast<T>((std::numeric_limits<boost::long_long_type>::min)()));
check_within_half(static_cast<T>((std::numeric_limits<boost::long_long_type>::min)()), j);
j = lltrunc(static_cast<T>((std::numeric_limits<boost::long_long_type>::max)()));
check_trunc_result(static_cast<T>((std::numeric_limits<boost::long_long_type>::max)()), j);
j = lltrunc(static_cast<T>((std::numeric_limits<boost::long_long_type>::min)()));
check_trunc_result(static_cast<T>((std::numeric_limits<boost::long_long_type>::min)()), j);
}
#endif
test_round_number<T>(arg);
}
//
// Finish off by testing the error handlers:
Expand Down Expand Up @@ -360,10 +371,66 @@ BOOST_AUTO_TEST_CASE( test_main )
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::cout;
#endif

}

// test rounding of direct predecessor/successor of 0.5/-0.5 for float and double
test_round_number(-0.4999999701976776123046875f);
BOOST_CHECK_EQUAL(boost::math::round(-0.4999999701976776123046875f), 0.0f);

test_round_number(0.4999999701976776123046875f);
BOOST_CHECK_EQUAL(boost::math::round(0.4999999701976776123046875f), 0.0f);

BOOST_CHECK_EQUAL(boost::math::round(-0.499999999999999944488848768742172978818416595458984375), 0.0);
test_round_number(-0.499999999999999944488848768742172978818416595458984375);

BOOST_CHECK_EQUAL(boost::math::round(0.499999999999999944488848768742172978818416595458984375), 0.0);
test_round_number(0.499999999999999944488848768742172978818416595458984375);

// test rounding of integer numbers on the edge of the float/double mantissa width
BOOST_CHECK_EQUAL(boost::math::round(-16777215.0f), -16777215.0f);
test_round_number(-16777215.0f);

BOOST_CHECK_EQUAL(boost::math::round(-16777213.0f), -16777213.0f);
test_round_number(-16777213.0f);

BOOST_CHECK_EQUAL(boost::math::round(-8388611.0f), -8388611.0f);
test_round_number(-8388611.0f);

BOOST_CHECK_EQUAL(boost::math::round(-8388609.0f), -8388609.0f);
test_round_number(-8388609.0f);

BOOST_CHECK_EQUAL(boost::math::round(8388609.0f), 8388609.0f);
test_round_number(8388609.0f);

BOOST_CHECK_EQUAL(boost::math::round(8388611.0f), 8388611.0f);
test_round_number(8388611.0f);

BOOST_CHECK_EQUAL(boost::math::round(16777213.0f), 16777213.0f);
test_round_number(16777213.0f);

BOOST_CHECK_EQUAL(boost::math::round(16777215.0f), 16777215.0f);
test_round_number(16777215.0f);

BOOST_CHECK_EQUAL(boost::math::round(-9007199254740993.0), -9007199254740993.0);
test_round_number(-9007199254740993.0);

BOOST_CHECK_EQUAL(boost::math::round(-9007199254740991.0), -9007199254740991.0);
test_round_number(-9007199254740991.0);

BOOST_CHECK_EQUAL(boost::math::round(-4503599627370499.0), -4503599627370499.0);
test_round_number(-4503599627370499.0);

BOOST_CHECK_EQUAL(boost::math::round(-4503599627370497.0), -4503599627370497.0);
test_round_number(-4503599627370497.0);

BOOST_CHECK_EQUAL(boost::math::round(4503599627370497.0), 4503599627370497.0);
test_round_number(4503599627370497.0);

BOOST_CHECK_EQUAL(boost::math::round(4503599627370499.0), 4503599627370499.0);
test_round_number(4503599627370499.0);

BOOST_CHECK_EQUAL(boost::math::round(9007199254740991.0), 9007199254740991.0);
test_round_number(9007199254740991.0);

BOOST_CHECK_EQUAL(boost::math::round(9007199254740993.0), 9007199254740993.0);
test_round_number(9007199254740993.0);
}