| // (C) Copyright John Maddock 2008. |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. (See accompanying file |
| // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| #ifndef BOOST_MATH_SPECIAL_NEXT_HPP |
| #define BOOST_MATH_SPECIAL_NEXT_HPP |
| |
| #ifdef _MSC_VER |
| #pragma once |
| #endif |
| |
| #include <boost/math/policies/error_handling.hpp> |
| #include <boost/math/special_functions/fpclassify.hpp> |
| #include <boost/math/special_functions/sign.hpp> |
| #include <boost/math/special_functions/trunc.hpp> |
| |
| #ifdef BOOST_MSVC |
| #include <float.h> |
| #endif |
| |
| namespace boost{ namespace math{ |
| |
| namespace detail{ |
| |
| template <class T> |
| inline T get_smallest_value(mpl::true_ const&) |
| { |
| return std::numeric_limits<T>::denorm_min(); |
| } |
| |
| template <class T> |
| inline T get_smallest_value(mpl::false_ const&) |
| { |
| return tools::min_value<T>(); |
| } |
| |
| template <class T> |
| inline T get_smallest_value() |
| { |
| #if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310) |
| return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>()); |
| #else |
| return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>()); |
| #endif |
| } |
| |
| } |
| |
| template <class T, class Policy> |
| T float_next(const T& val, const Policy& pol) |
| { |
| BOOST_MATH_STD_USING |
| int expon; |
| static const char* function = "float_next<%1%>(%1%)"; |
| |
| if(!(boost::math::isfinite)(val)) |
| 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(-0.5f == frexp(val, &expon)) |
| --expon; // reduce exponent when val is a power of two, and negative. |
| T diff = ldexp(T(1), expon - tools::digits<T>()); |
| if(diff == 0) |
| diff = detail::get_smallest_value<T>(); |
| return val + diff; |
| } |
| |
| #ifdef BOOST_MSVC |
| template <class Policy> |
| inline double float_next(const double& val, const Policy& pol) |
| { |
| static const char* function = "float_next<%1%>(%1%)"; |
| |
| if(!(boost::math::isfinite)(val)) |
| return policies::raise_domain_error<double>( |
| function, |
| "Argument must be finite, but got %1%", val, pol); |
| |
| if(val >= tools::max_value<double>()) |
| return policies::raise_overflow_error<double>(function, 0, pol); |
| |
| return ::_nextafter(val, tools::max_value<double>()); |
| } |
| #endif |
| |
| template <class T> |
| inline T float_next(const T& val) |
| { |
| return float_next(val, policies::policy<>()); |
| } |
| |
| template <class T, class Policy> |
| T float_prior(const T& val, const Policy& pol) |
| { |
| BOOST_MATH_STD_USING |
| int expon; |
| static const char* function = "float_prior<%1%>(%1%)"; |
| |
| if(!(boost::math::isfinite)(val)) |
| 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>(); |
| |
| T remain = frexp(val, &expon); |
| if(remain == 0.5) |
| --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; |
| } |
| |
| #ifdef BOOST_MSVC |
| template <class Policy> |
| inline double float_prior(const double& val, const Policy& pol) |
| { |
| static const char* function = "float_prior<%1%>(%1%)"; |
| |
| if(!(boost::math::isfinite)(val)) |
| return policies::raise_domain_error<double>( |
| function, |
| "Argument must be finite, but got %1%", val, pol); |
| |
| if(val <= -tools::max_value<double>()) |
| return -policies::raise_overflow_error<double>(function, 0, pol); |
| |
| return ::_nextafter(val, -tools::max_value<double>()); |
| } |
| #endif |
| |
| template <class T> |
| inline T float_prior(const T& val) |
| { |
| return float_prior(val, policies::policy<>()); |
| } |
| |
| template <class T, class Policy> |
| inline T nextafter(const T& val, const T& direction, const Policy& pol) |
| { |
| return val < direction ? boost::math::float_next(val, pol) : val == direction ? val : boost::math::float_prior(val, pol); |
| } |
| |
| template <class T> |
| inline T nextafter(const T& val, const T& direction) |
| { |
| return nextafter(val, direction, policies::policy<>()); |
| } |
| |
| template <class T, class Policy> |
| T float_distance(const T& a, const T& b, const Policy& pol) |
| { |
| 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); |
| if(a == b) |
| return 0; |
| if(a == 0) |
| return 1 + fabs(float_distance(static_cast<T>(boost::math::sign(b) * detail::get_smallest_value<T>()), b, pol)); |
| if(b == 0) |
| return 1 + fabs(float_distance(static_cast<T>(boost::math::sign(a) * detail::get_smallest_value<T>()), a, pol)); |
| if(boost::math::sign(a) != boost::math::sign(b)) |
| return 2 + fabs(float_distance(static_cast<T>(boost::math::sign(b) * detail::get_smallest_value<T>()), b, pol)) |
| + fabs(float_distance(static_cast<T>(boost::math::sign(a) * 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)); |
| |
| BOOST_ASSERT(a >= 0); |
| BOOST_ASSERT(b >= a); |
| |
| int 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: |
| // |
| frexp(((boost::math::fpclassify)(a) == FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon); |
| T upper = ldexp(T(1), expon); |
| T result = 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); |
| } |
| // |
| // Use compensated double-double addition to avoid rounding |
| // errors in the subtraction: |
| // |
| T mb = -(std::min)(upper, b); |
| T x = a + mb; |
| T z = x - a; |
| T y = (a - (x - z)) + (mb - z); |
| if(x < 0) |
| { |
| x = -x; |
| y = -y; |
| } |
| result += ldexp(x, expon) + ldexp(y, expon); |
| // |
| // Result must be an integer: |
| // |
| BOOST_ASSERT(result == floor(result)); |
| return result; |
| } |
| |
| template <class T> |
| T float_distance(const T& a, const T& b) |
| { |
| return boost::math::float_distance(a, b, policies::policy<>()); |
| } |
| |
| template <class T, class Policy> |
| T float_advance(T val, int distance, const Policy& pol) |
| { |
| // |
| // Error handling: |
| // |
| static const char* function = "float_advance<%1%>(%1%, int)"; |
| if(!(boost::math::isfinite)(val)) |
| 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); |
| BOOST_MATH_STD_USING |
| int expon; |
| frexp(val, &expon); |
| T limit = ldexp((distance < 0 ? T(0.5f) : T(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 /= 2; |
| expon--; |
| } |
| else |
| { |
| limit *= 2; |
| expon++; |
| } |
| limit_distance = float_distance(val, limit); |
| } |
| if((0.5f == frexp(val, &expon)) && (distance < 0)) |
| --expon; |
| T diff = 0; |
| if(val != 0) |
| diff = distance * ldexp(T(1), expon - tools::digits<T>()); |
| if(diff == 0) |
| diff = distance * detail::get_smallest_value<T>(); |
| return val += diff; |
| } |
| |
| template <class T> |
| inline T float_advance(const T& val, int distance) |
| { |
| return boost::math::float_advance(val, distance, policies::policy<>()); |
| } |
| |
| }} // namespaces |
| |
| #endif // BOOST_MATH_SPECIAL_NEXT_HPP |
| |