346 changes: 336 additions & 10 deletions include/boost/math/special_functions/next.hpp
View file
Expand Up @@ -29,6 +29,28 @@ namespace boost{ namespace math{

namespace detail{

template <class T>
struct has_hidden_guard_digits
: public mpl::bool_
<std::numeric_limits<T>::is_specialized
&& (std::numeric_limits<T>::radix == 10)
&& (std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)>
{};

template <class T>
inline const T& normalize_value(const T& val, const mpl::false_&) { return val; }
template <class T>
inline T normalize_value(const T& val, const mpl::true_&)
{
BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);

boost::intmax_t shift = std::numeric_limits<T>::digits - ilogb(val) - 1;
T result = scalbn(val, shift);
result = round(result);
return scalbn(result, -shift);
}

template <class T>
inline T get_smallest_value(mpl::true_ const&)
{
Expand Down Expand Up @@ -93,19 +115,33 @@ struct min_shift_initializer
template <class T>
const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;

template <class T>
inline T calc_min_shifted(const mpl::true_&)
{
BOOST_MATH_STD_USING
return ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
}
template <class T>
inline T calc_min_shifted(const mpl::false_&)
{
BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);

return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);
}


template <class T>
inline T get_min_shift_value()
{
BOOST_MATH_STD_USING
static const T val = ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
static const T val = calc_min_shifted<T>(mpl::bool_<!std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());
min_shift_initializer<T>::force_instantiate();

return val;
}

template <class T, class Policy>
T float_next_imp(const T& val, const Policy& pol)
T float_next_imp(const T& val, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
int expon;
Expand Down Expand Up @@ -145,14 +181,62 @@ T float_next_imp(const T& val, const Policy& pol)
diff = detail::get_smallest_value<T>();
return val + diff;
} // float_next_imp
//
// Special version for some base other than 2:
//
template <class T, class Policy>
T float_next_imp(const T& val, const mpl::false_&, const Policy& pol)
{
BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);

BOOST_MATH_STD_USING
boost::intmax_t expon;
static const char* function = "float_next<%1%>(%1%)";

int fpclass = (boost::math::fpclassify)(val);

if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
{
if(val < 0)
return -tools::max_value<T>();
return policies::raise_domain_error<T>(
function,
"Argument must be finite, but got %1%", val, pol);
}

if(val >= tools::max_value<T>())
return policies::raise_overflow_error<T>(function, 0, pol);

if(val == 0)
return detail::get_smallest_value<T>();

if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
{
//
// Special case: if the value of the least significant bit is a denorm, and the result
// would not be a denorm, then shift the input, increment, and shift back.
// This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
//
return scalbn(float_next(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
}

expon = 1 + ilogb(val);
if(-1 == scalbn(val, -expon) * std::numeric_limits<T>::radix)
--expon; // reduce exponent when val is a power of base, and negative.
T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
if(diff == 0)
diff = detail::get_smallest_value<T>();
return val + diff;
} // float_next_imp

} // namespace detail

template <class T, class Policy>
inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
return detail::float_next_imp(static_cast<result_type>(val), pol);
return detail::float_next_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<T>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}

#if 0 //def BOOST_MSVC
Expand Down Expand Up @@ -187,7 +271,7 @@ inline typename tools::promote_args<T>::type float_next(const T& val)
namespace detail{

template <class T, class Policy>
T float_prior_imp(const T& val, const Policy& pol)
T float_prior_imp(const T& val, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
int expon;
Expand Down Expand Up @@ -228,14 +312,63 @@ T float_prior_imp(const T& val, const Policy& pol)
diff = detail::get_smallest_value<T>();
return val - diff;
} // float_prior_imp
//
// Special version for bases other than 2:
//
template <class T, class Policy>
T float_prior_imp(const T& val, const mpl::false_&, const Policy& pol)
{
BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);

BOOST_MATH_STD_USING
boost::intmax_t expon;
static const char* function = "float_prior<%1%>(%1%)";

int fpclass = (boost::math::fpclassify)(val);

if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
{
if(val > 0)
return tools::max_value<T>();
return policies::raise_domain_error<T>(
function,
"Argument must be finite, but got %1%", val, pol);
}

if(val <= -tools::max_value<T>())
return -policies::raise_overflow_error<T>(function, 0, pol);

if(val == 0)
return -detail::get_smallest_value<T>();

if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
{
//
// Special case: if the value of the least significant bit is a denorm, and the result
// would not be a denorm, then shift the input, increment, and shift back.
// This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
//
return scalbn(float_prior(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
}

expon = 1 + ilogb(val);
T remain = scalbn(val, -expon);
if(remain * std::numeric_limits<T>::radix == 1)
--expon; // when val is a power of two we must reduce the exponent
T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
if(diff == 0)
diff = detail::get_smallest_value<T>();
return val - diff;
} // float_prior_imp

} // namespace detail

template <class T, class Policy>
inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
return detail::float_prior_imp(static_cast<result_type>(val), pol);
return detail::float_prior_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<T>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}

#if 0 //def BOOST_MSVC
Expand Down Expand Up @@ -283,7 +416,7 @@ inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U&
namespace detail{

template <class T, class Policy>
T float_distance_imp(const T& a, const T& b, const Policy& pol)
T float_distance_imp(const T& a, const T& b, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
//
Expand Down Expand Up @@ -384,14 +517,121 @@ T float_distance_imp(const T& a, const T& b, const Policy& pol)
BOOST_ASSERT(result == floor(result));
return result;
} // float_distance_imp
//
// Special versions for bases other than 2:
//
template <class T, class Policy>
T float_distance_imp(const T& a, const T& b, const mpl::false_&, const Policy& pol)
{
BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);

BOOST_MATH_STD_USING
//
// Error handling:
//
static const char* function = "float_distance<%1%>(%1%, %1%)";
if(!(boost::math::isfinite)(a))
return policies::raise_domain_error<T>(
function,
"Argument a must be finite, but got %1%", a, pol);
if(!(boost::math::isfinite)(b))
return policies::raise_domain_error<T>(
function,
"Argument b must be finite, but got %1%", b, pol);
//
// Special cases:
//
if(a > b)
return -float_distance(b, a, pol);
if(a == b)
return T(0);
if(a == 0)
return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
if(b == 0)
return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
if(boost::math::sign(a) != boost::math::sign(b))
return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
+ fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
//
// By the time we get here, both a and b must have the same sign, we want
// b > a and both postive for the following logic:
//
if(a < 0)
return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);

BOOST_ASSERT(a >= 0);
BOOST_ASSERT(b >= a);

boost::intmax_t expon;
//
// Note that if a is a denorm then the usual formula fails
// because we actually have fewer than tools::digits<T>()
// significant bits in the representation:
//
expon = 1 + ilogb(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a);
T upper = scalbn(T(1), expon);
T result = T(0);
//
// If b is greater than upper, then we *must* split the calculation
// as the size of the ULP changes with each order of magnitude change:
//
if(b > upper)
{
boost::intmax_t expon2 = 1 + ilogb(b);
T upper2 = scalbn(T(1), expon2 - 1);
result = float_distance(upper2, b);
result += (expon2 - expon - 1) * scalbn(T(1), std::numeric_limits<T>::digits - 1);
}
//
// Use compensated double-double addition to avoid rounding
// errors in the subtraction:
//
expon = std::numeric_limits<T>::digits - expon;
T mb, x, y, z;
if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
{
//
// Special case - either one end of the range is a denormal, or else the difference is.
// The regular code will fail if we're using the SSE2 registers on Intel and either
// the FTZ or DAZ flags are set.
//
T a2 = scalbn(a, std::numeric_limits<T>::digits);
T b2 = scalbn(b, std::numeric_limits<T>::digits);
mb = -(std::min)(T(scalbn(upper, std::numeric_limits<T>::digits)), b2);
x = a2 + mb;
z = x - a2;
y = (a2 - (x - z)) + (mb - z);

expon -= std::numeric_limits<T>::digits;
}
else
{
mb = -(std::min)(upper, b);
x = a + mb;
z = x - a;
y = (a - (x - z)) + (mb - z);
}
if(x < 0)
{
x = -x;
y = -y;
}
result += scalbn(x, expon) + scalbn(y, expon);
//
// Result must be an integer:
//
BOOST_ASSERT(result == floor(result));
return result;
} // float_distance_imp

} // namespace detail

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
{
typedef typename tools::promote_args<T, U>::type result_type;
return detail::float_distance_imp(static_cast<result_type>(a), static_cast<result_type>(b), pol);
return detail::float_distance_imp(detail::normalize_value(static_cast<result_type>(a), typename detail::has_hidden_guard_digits<T>::type()), detail::normalize_value(static_cast<result_type>(b), typename detail::has_hidden_guard_digits<T>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}

template <class T, class U>
Expand All @@ -403,7 +643,7 @@ typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
namespace detail{

template <class T, class Policy>
T float_advance_imp(T val, int distance, const Policy& pol)
T float_advance_imp(T val, int distance, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
//
Expand Down Expand Up @@ -482,14 +722,100 @@ T float_advance_imp(T val, int distance, const Policy& pol)
diff = distance * detail::get_smallest_value<T>();
return val += diff;
} // float_advance_imp
//
// Special version for bases other than 2:
//
template <class T, class Policy>
T float_advance_imp(T val, int distance, const mpl::false_&, const Policy& pol)
{
BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);

BOOST_MATH_STD_USING
//
// Error handling:
//
static const char* function = "float_advance<%1%>(%1%, int)";

int fpclass = (boost::math::fpclassify)(val);

if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
return policies::raise_domain_error<T>(
function,
"Argument val must be finite, but got %1%", val, pol);

if(val < 0)
return -float_advance(-val, -distance, pol);
if(distance == 0)
return val;
if(distance == 1)
return float_next(val, pol);
if(distance == -1)
return float_prior(val, pol);

if(fabs(val) < detail::get_min_shift_value<T>())
{
//
// Special case: if the value of the least significant bit is a denorm,
// implement in terms of float_next/float_prior.
// This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
//
if(distance > 0)
{
do{ val = float_next(val, pol); } while(--distance);
}
else
{
do{ val = float_prior(val, pol); } while(++distance);
}
return val;
}

boost::intmax_t expon = 1 + ilogb(val);
T limit = scalbn(T(1), distance < 0 ? expon - 1 : expon);
if(val <= tools::min_value<T>())
{
limit = sign(T(distance)) * tools::min_value<T>();
}
T limit_distance = float_distance(val, limit);
while(fabs(limit_distance) < abs(distance))
{
distance -= itrunc(limit_distance);
val = limit;
if(distance < 0)
{
limit /= std::numeric_limits<T>::radix;
expon--;
}
else
{
limit *= std::numeric_limits<T>::radix;
expon++;
}
limit_distance = float_distance(val, limit);
if(distance && (limit_distance == 0))
{
return policies::raise_evaluation_error<T>(function, "Internal logic failed while trying to increment floating point value %1%: most likely your FPU is in non-IEEE conforming mode.", val, pol);
}
}
/*expon = 1 + ilogb(val);
if((1 == scalbn(val, 1 + expon)) && (distance < 0))
--expon;*/
T diff = 0;
if(val != 0)
diff = distance * scalbn(T(1), expon - std::numeric_limits<T>::digits);
if(diff == 0)
diff = distance * detail::get_smallest_value<T>();
return val += diff;
} // float_advance_imp

} // namespace detail

template <class T, class Policy>
inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
return detail::float_advance_imp(static_cast<result_type>(val), distance, pol);
return detail::float_advance_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<T>::type()), distance, mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}

template <class T>
Expand Down
38 changes: 36 additions & 2 deletions include/boost/math/special_functions/ulp.hpp
View file
Expand Up @@ -18,7 +18,7 @@
namespace boost{ namespace math{ namespace detail{

template <class T, class Policy>
T ulp_imp(const T& val, const Policy& pol)
T ulp_imp(const T& val, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
int expon;
Expand Down Expand Up @@ -48,14 +48,48 @@ T ulp_imp(const T& val, const Policy& pol)
diff = detail::get_smallest_value<T>();
return diff;
}
// non-binary version:
template <class T, class Policy>
T ulp_imp(const T& val, const mpl::false_&, const Policy& pol)
{
BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
BOOST_MATH_STD_USING
int expon;
static const char* function = "ulp<%1%>(%1%)";

int fpclass = (boost::math::fpclassify)(val);

if(fpclass == (int)FP_NAN)
{
return policies::raise_domain_error<T>(
function,
"Argument must be finite, but got %1%", val, pol);
}
else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
{
return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, 0, pol);
}
else if(fpclass == FP_ZERO)
return detail::get_smallest_value<T>();
//
// This code is almost the same as that for float_next, except for negative integers,
// where we preserve the relation ulp(x) == ulp(-x) as does Java:
//
expon = 1 + ilogb(fabs(val));
T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
if(diff == 0)
diff = detail::get_smallest_value<T>();
return diff;
}

}

template <class T, class Policy>
inline typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
return detail::ulp_imp(static_cast<result_type>(val), pol);
return detail::ulp_imp(static_cast<result_type>(val), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}

template <class T>
Expand Down
2 changes: 2 additions & 0 deletions test/Jamfile.v2
View file
Expand Up @@ -598,6 +598,7 @@ run test_negative_binomial.cpp ../../test/build//boost_unit_test_framework
: test_negative_binomial_real_concept ;
# Needs to run in release mode, as it's rather slow:
run test_next.cpp pch ../../test/build//boost_unit_test_framework : : : release ;
run test_next_decimal.cpp pch ../../test/build//boost_unit_test_framework : : : release ;
run test_nc_chi_squared.cpp pch ../../test/build//boost_unit_test_framework
: # command line
: # input files
Expand Down Expand Up @@ -1173,6 +1174,7 @@ legendre_stieltjes_test
test_ldouble_simple
test_minima
test_next
test_next_decimal
test_rationals
test_real_concept
test_remez
Expand Down