diff options
Diffstat (limited to 'boost/math/special_functions/next.hpp')
| -rw-r--r-- | boost/math/special_functions/next.hpp | 384 |
1 files changed, 371 insertions, 13 deletions
diff --git a/boost/math/special_functions/next.hpp b/boost/math/special_functions/next.hpp index 3921b70ce0..f23ddaffe3 100644 --- a/boost/math/special_functions/next.hpp +++ b/boost/math/special_functions/next.hpp @@ -27,9 +27,59 @@ namespace boost{ namespace math{ + namespace concepts { + + struct real_concept; + struct std_real_concept; + + } + namespace detail{ template <class T> +struct has_hidden_guard_digits; +template <> +struct has_hidden_guard_digits<float> : public mpl::false_ {}; +template <> +struct has_hidden_guard_digits<double> : public mpl::false_ {}; +template <> +struct has_hidden_guard_digits<long double> : public mpl::false_ {}; +#ifdef BOOST_HAS_FLOAT128 +template <> +struct has_hidden_guard_digits<__float128> : public mpl::false_ {}; +#endif +template <> +struct has_hidden_guard_digits<boost::math::concepts::real_concept> : public mpl::false_ {}; +template <> +struct has_hidden_guard_digits<boost::math::concepts::std_real_concept> : public mpl::false_ {}; + +template <class T, bool b> +struct has_hidden_guard_digits_10 : public mpl::false_ {}; +template <class T> +struct has_hidden_guard_digits_10<T, true> : public mpl::bool_<(std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {}; + +template <class T> +struct has_hidden_guard_digits + : public has_hidden_guard_digits_10<T, + std::numeric_limits<T>::is_specialized + && (std::numeric_limits<T>::radix == 10) > +{}; + +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&) { // @@ -93,19 +143,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; @@ -145,6 +209,54 @@ 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 @@ -152,7 +264,7 @@ 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<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); } #if 0 //def BOOST_MSVC @@ -187,7 +299,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; @@ -221,13 +333,62 @@ T float_prior_imp(const T& val, const Policy& pol) } T remain = frexp(val, &expon); - if(remain == 0.5) + if(remain == 0.5f) --expon; // when val is a power of two we must reduce the exponent T diff = ldexp(T(1), expon - tools::digits<T>()); if(diff == 0) 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 @@ -235,7 +396,7 @@ 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<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); } #if 0 //def BOOST_MSVC @@ -283,7 +444,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 // @@ -331,19 +492,23 @@ T float_distance_imp(const T& a, const T& b, const Policy& pol) frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon); T upper = ldexp(T(1), expon); T result = T(0); - expon = tools::digits<T>() - expon; // // 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) { - result = float_distance(upper, b); + int expon2; + frexp(b, &expon2); + T upper2 = ldexp(T(0.5), expon2); + result = float_distance(upper2, b); + result += (expon2 - expon - 1) * ldexp(T(1), tools::digits<T>() - 1); } // // Use compensated double-double addition to avoid rounding // errors in the subtraction: // + expon = tools::digits<T>() - expon; T mb, x, y, z; if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>())) { @@ -380,6 +545,113 @@ 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 @@ -387,7 +659,7 @@ 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<result_type>::type()), detail::normalize_value(static_cast<result_type>(b), typename detail::has_hidden_guard_digits<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); } template <class T, class U> @@ -399,7 +671,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 // @@ -478,6 +750,92 @@ 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 @@ -485,7 +843,7 @@ 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<result_type>::type()), distance, mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); } template <class T> |