math/special_functions/next.hpp

0001 //  (C) Copyright John Maddock 2008.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SPECIAL_NEXT_HPP
0007 #define BOOST_MATH_SPECIAL_NEXT_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/config.hpp>
0014
0015 // TODO(mborland): Need to remove recurrsion from these algos
0016 #ifndef BOOST_MATH_HAS_NVRTC
0017
0018 #include <boost/math/special_functions/math_fwd.hpp>
0019 #include <boost/math/policies/error_handling.hpp>
0020 #include <boost/math/special_functions/fpclassify.hpp>
0021 #include <boost/math/special_functions/sign.hpp>
0022 #include <boost/math/special_functions/trunc.hpp>
0023 #include <boost/math/tools/traits.hpp>
0024 #include <type_traits>
0025 #include <cfloat>
0026
0027
0028 #if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3)))
0029 #if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__)
0030 #include "xmmintrin.h"
0031 #define BOOST_MATH_CHECK_SSE2
0032 #endif
0033 #endif
0034
0035 namespace boost{ namespace math{
0036
0037    namespace concepts {
0038
0039       class real_concept;
0040       class std_real_concept;
0041
0042    }
0043
0044 namespace detail{
0045
0046 template <class T>
0047 struct has_hidden_guard_digits;
0048 template <>
0049 struct has_hidden_guard_digits<float> : public std::false_type {};
0050 template <>
0051 struct has_hidden_guard_digits<double> : public std::false_type {};
0052 template <>
0053 struct has_hidden_guard_digits<long double> : public std::false_type {};
0054 #ifdef BOOST_HAS_FLOAT128
0055 template <>
0056 struct has_hidden_guard_digits<__float128> : public std::false_type {};
0057 #endif
0058 template <>
0059 struct has_hidden_guard_digits<boost::math::concepts::real_concept> : public std::false_type {};
0060 template <>
0061 struct has_hidden_guard_digits<boost::math::concepts::std_real_concept> : public std::false_type {};
0062
0063 template <class T, bool b>
0064 struct has_hidden_guard_digits_10 : public std::false_type {};
0065 template <class T>
0066 struct has_hidden_guard_digits_10<T, true> : public std::integral_constant<bool, (std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {};
0067
0068 template <class T>
0069 struct has_hidden_guard_digits
0070    : public has_hidden_guard_digits_10<T,
0071    std::numeric_limits<T>::is_specialized
0072    && (std::numeric_limits<T>::radix == 10) >
0073 {};
0074
0075 template <class T>
0076 inline const T& normalize_value(const T& val, const std::false_type&) { return val; }
0077 template <class T>
0078 inline T normalize_value(const T& val, const std::true_type&)
0079 {
0080    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0081    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0082
0083    std::intmax_t shift = (std::intmax_t)std::numeric_limits<T>::digits - (std::intmax_t)ilogb(val) - 1;
0084    T result = scalbn(val, shift);
0085    result = round(result);
0086    return scalbn(result, -shift);
0087 }
0088
0089 template <class T>
0090 inline T get_smallest_value(std::true_type const&) {
0091    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0092    //
0093    // numeric_limits lies about denorms being present - particularly
0094    // when this can be turned on or off at runtime, as is the case
0095    // when using the SSE2 registers in DAZ or FTZ mode.
0096    //
0097    static const T m = std::numeric_limits<T>::denorm_min();
0098 #ifdef BOOST_MATH_CHECK_SSE2
0099    return (_mm_getcsr() & (_MM_FLUSH_ZERO_ON | 0x40)) ? tools::min_value<T>() : m;
0100 #else
0101    return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m;
0102 #endif
0103 }
0104
0105 template <class T>
0106 inline T get_smallest_value(std::false_type const&)
0107 {
0108    return tools::min_value<T>();
0109 }
0110
0111 template <class T>
0112 inline T get_smallest_value()
0113 {
0114    return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized>());
0115 }
0116
0117 template <class T>
0118 inline bool has_denorm_now() {
0119    return get_smallest_value<T>() < tools::min_value<T>();
0120 }
0121
0122 //
0123 // Returns the smallest value that won't generate denorms when
0124 // we calculate the value of the least-significant-bit:
0125 //
0126 template <class T>
0127 T get_min_shift_value();
0128
0129 template <class T>
0130 inline T calc_min_shifted(const std::true_type&)
0131 {
0132    BOOST_MATH_STD_USING
0133    return ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
0134 }
0135 template <class T>
0136 inline T calc_min_shifted(const std::false_type&)
0137 {
0138    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0139    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0140
0141    return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);
0142 }
0143
0144
0145 template <class T>
0146 inline T get_min_shift_value()
0147 {
0148    static const T val = calc_min_shifted<T>(std::integral_constant<bool, !std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());
0149    return val;
0150 }
0151
0152 template <class T, bool b = boost::math::tools::detail::has_backend_type<T>::value>
0153 struct exponent_type
0154 {
0155    typedef int type;
0156 };
0157
0158 template <class T>
0159 struct exponent_type<T, true>
0160 {
0161    typedef typename T::backend_type::exponent_type type;
0162 };
0163
0164 template <class T, class Policy>
0165 T float_next_imp(const T& val, const std::true_type&, const Policy& pol)
0166 {
0167    typedef typename exponent_type<T>::type exponent_type;
0168
0169    BOOST_MATH_STD_USING
0170    exponent_type expon;
0171    static const char* function = "float_next<%1%>(%1%)";
0172
0173    int fpclass = (boost::math::fpclassify)(val);
0174
0175    if (fpclass == (int)FP_INFINITE)
0176    {
0177       if (val < 0)
0178          return -tools::max_value<T>();
0179       return val;  // +INF
0180    }
0181    else if (fpclass == (int)FP_NAN)
0182    {
0183       return policies::raise_domain_error<T>(
0184          function,
0185          "Argument must be finite, but got %1%", val, pol);
0186    }
0187
0188    if(val >= tools::max_value<T>())
0189       return policies::raise_overflow_error<T>(function, nullptr, pol);
0190
0191    if(val == 0)
0192       return detail::get_smallest_value<T>();
0193
0194    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
0195    {
0196       //
0197       // Special case: if the value of the least significant bit is a denorm, and the result
0198       // would not be a denorm, then shift the input, increment, and shift back.
0199       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0200       //
0201       return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
0202    }
0203
0204    if(-0.5f == frexp(val, &expon))
0205       --expon; // reduce exponent when val is a power of two, and negative.
0206    T diff = ldexp(T(1), expon - tools::digits<T>());
0207    if(diff == 0)
0208       diff = detail::get_smallest_value<T>();
0209    return val + diff;
0210 } // float_next_imp
0211 //
0212 // Special version for some base other than 2:
0213 //
0214 template <class T, class Policy>
0215 T float_next_imp(const T& val, const std::false_type&, const Policy& pol)
0216 {
0217    typedef typename exponent_type<T>::type exponent_type;
0218
0219    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0220    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0221
0222    BOOST_MATH_STD_USING
0223    exponent_type expon;
0224    static const char* function = "float_next<%1%>(%1%)";
0225
0226    int fpclass = (boost::math::fpclassify)(val);
0227
0228    if (fpclass == (int)FP_INFINITE)
0229    {
0230       if (val < 0)
0231          return -tools::max_value<T>();
0232       return val;  // +INF
0233    }
0234    else if (fpclass == (int)FP_NAN)
0235    {
0236       return policies::raise_domain_error<T>(
0237          function,
0238          "Argument must be finite, but got %1%", val, pol);
0239    }
0240
0241    if(val >= tools::max_value<T>())
0242       return policies::raise_overflow_error<T>(function, nullptr, pol);
0243
0244    if(val == 0)
0245       return detail::get_smallest_value<T>();
0246
0247    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
0248    {
0249       //
0250       // Special case: if the value of the least significant bit is a denorm, and the result
0251       // would not be a denorm, then shift the input, increment, and shift back.
0252       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0253       //
0254       return scalbn(float_next(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
0255    }
0256
0257    expon = 1 + ilogb(val);
0258    if(-1 == scalbn(val, -expon) * std::numeric_limits<T>::radix)
0259       --expon; // reduce exponent when val is a power of base, and negative.
0260    T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
0261    if(diff == 0)
0262       diff = detail::get_smallest_value<T>();
0263    return val + diff;
0264 } // float_next_imp
0265
0266 } // namespace detail
0267
0268 template <class T, class Policy>
0269 inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
0270 {
0271    typedef typename tools::promote_args<T>::type result_type;
0272    return detail::float_next_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
0273 }
0274
0275 #if 0 //def BOOST_MSVC
0276 //
0277 // We used to use ::_nextafter here, but doing so fails when using
0278 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
0279 // - albeit slower - code instead as at least that gives the correct answer.
0280 //
0281 template <class Policy>
0282 inline double float_next(const double& val, const Policy& pol)
0283 {
0284    static const char* function = "float_next<%1%>(%1%)";
0285
0286    if(!(boost::math::isfinite)(val) && (val > 0))
0287       return policies::raise_domain_error<double>(
0288          function,
0289          "Argument must be finite, but got %1%", val, pol);
0290
0291    if(val >= tools::max_value<double>())
0292       return policies::raise_overflow_error<double>(function, nullptr, pol);
0293
0294    return ::_nextafter(val, tools::max_value<double>());
0295 }
0296 #endif
0297
0298 template <class T>
0299 inline typename tools::promote_args<T>::type float_next(const T& val)
0300 {
0301    return float_next(val, policies::policy<>());
0302 }
0303
0304 namespace detail{
0305
0306 template <class T, class Policy>
0307 T float_prior_imp(const T& val, const std::true_type&, const Policy& pol)
0308 {
0309    typedef typename exponent_type<T>::type exponent_type;
0310
0311    BOOST_MATH_STD_USING
0312    exponent_type expon;
0313    static const char* function = "float_prior<%1%>(%1%)";
0314
0315    int fpclass = (boost::math::fpclassify)(val);
0316
0317    if (fpclass == (int)FP_INFINITE)
0318    {
0319       if (val > 0)
0320          return tools::max_value<T>();
0321       return val; // -INF
0322    }
0323    else if (fpclass == (int)FP_NAN)
0324    {
0325       return policies::raise_domain_error<T>(
0326          function,
0327          "Argument must be finite, but got %1%", val, pol);
0328    }
0329
0330    if(val <= -tools::max_value<T>())
0331       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0332
0333    if(val == 0)
0334       return -detail::get_smallest_value<T>();
0335
0336    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
0337    {
0338       //
0339       // Special case: if the value of the least significant bit is a denorm, and the result
0340       // would not be a denorm, then shift the input, increment, and shift back.
0341       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0342       //
0343       return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
0344    }
0345
0346    T remain = frexp(val, &expon);
0347    if(remain == 0.5f)
0348       --expon; // when val is a power of two we must reduce the exponent
0349    T diff = ldexp(T(1), expon - tools::digits<T>());
0350    if(diff == 0)
0351       diff = detail::get_smallest_value<T>();
0352    return val - diff;
0353 } // float_prior_imp
0354 //
0355 // Special version for bases other than 2:
0356 //
0357 template <class T, class Policy>
0358 T float_prior_imp(const T& val, const std::false_type&, const Policy& pol)
0359 {
0360    typedef typename exponent_type<T>::type exponent_type;
0361
0362    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0363    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0364
0365    BOOST_MATH_STD_USING
0366    exponent_type expon;
0367    static const char* function = "float_prior<%1%>(%1%)";
0368
0369    int fpclass = (boost::math::fpclassify)(val);
0370
0371    if (fpclass == (int)FP_INFINITE)
0372    {
0373       if (val > 0)
0374          return tools::max_value<T>();
0375       return val; // -INF
0376    }
0377    else if (fpclass == (int)FP_NAN)
0378    {
0379       return policies::raise_domain_error<T>(
0380          function,
0381          "Argument must be finite, but got %1%", val, pol);
0382    }
0383
0384    if(val <= -tools::max_value<T>())
0385       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0386
0387    if(val == 0)
0388       return -detail::get_smallest_value<T>();
0389
0390    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
0391    {
0392       //
0393       // Special case: if the value of the least significant bit is a denorm, and the result
0394       // would not be a denorm, then shift the input, increment, and shift back.
0395       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0396       //
0397       return scalbn(float_prior(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
0398    }
0399
0400    expon = 1 + ilogb(val);
0401    T remain = scalbn(val, -expon);
0402    if(remain * std::numeric_limits<T>::radix == 1)
0403       --expon; // when val is a power of two we must reduce the exponent
0404    T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
0405    if(diff == 0)
0406       diff = detail::get_smallest_value<T>();
0407    return val - diff;
0408 } // float_prior_imp
0409
0410 } // namespace detail
0411
0412 template <class T, class Policy>
0413 inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
0414 {
0415    typedef typename tools::promote_args<T>::type result_type;
0416    return detail::float_prior_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
0417 }
0418
0419 #if 0 //def BOOST_MSVC
0420 //
0421 // We used to use ::_nextafter here, but doing so fails when using
0422 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
0423 // - albeit slower - code instead as at least that gives the correct answer.
0424 //
0425 template <class Policy>
0426 inline double float_prior(const double& val, const Policy& pol)
0427 {
0428    static const char* function = "float_prior<%1%>(%1%)";
0429
0430    if(!(boost::math::isfinite)(val) && (val < 0))
0431       return policies::raise_domain_error<double>(
0432          function,
0433          "Argument must be finite, but got %1%", val, pol);
0434
0435    if(val <= -tools::max_value<double>())
0436       return -policies::raise_overflow_error<double>(function, nullptr, pol);
0437
0438    return ::_nextafter(val, -tools::max_value<double>());
0439 }
0440 #endif
0441
0442 template <class T>
0443 inline typename tools::promote_args<T>::type float_prior(const T& val)
0444 {
0445    return float_prior(val, policies::policy<>());
0446 }
0447
0448 template <class T, class U, class Policy>
0449 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol)
0450 {
0451    typedef typename tools::promote_args<T, U>::type result_type;
0452    return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol);
0453 }
0454
0455 template <class T, class U>
0456 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
0457 {
0458    return nextafter(val, direction, policies::policy<>());
0459 }
0460
0461 namespace detail{
0462
0463 template <class T, class Policy>
0464 T float_distance_imp(const T& a, const T& b, const std::true_type&, const Policy& pol)
0465 {
0466    BOOST_MATH_STD_USING
0467    //
0468    // Error handling:
0469    //
0470    static const char* function = "float_distance<%1%>(%1%, %1%)";
0471    if(!(boost::math::isfinite)(a))
0472       return policies::raise_domain_error<T>(function, "Argument a must be finite, but got %1%", a, pol);
0473    if(!(boost::math::isfinite)(b))
0474       return policies::raise_domain_error<T>(function, "Argument b must be finite, but got %1%", b, pol);
0475    //
0476    // Special cases:
0477    //
0478    if(a > b)
0479       return -float_distance(b, a, pol);
0480    if(a == b)
0481       return T(0);
0482    if(a == 0)
0483       return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
0484    if(b == 0)
0485       return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
0486    if(boost::math::sign(a) != boost::math::sign(b))
0487       return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
0488          + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
0489    //
0490    // By the time we get here, both a and b must have the same sign, we want
0491    // b > a and both positive for the following logic:
0492    //
0493    if(a < 0)
0494       return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
0495
0496    BOOST_MATH_ASSERT(a >= 0);
0497    BOOST_MATH_ASSERT(b >= a);
0498
0499    int expon;
0500    //
0501    // Note that if a is a denorm then the usual formula fails
0502    // because we actually have fewer than tools::digits<T>()
0503    // significant bits in the representation:
0504    //
0505    (void)frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
0506    T upper = ldexp(T(1), expon);
0507    T result = T(0);
0508    //
0509    // If b is greater than upper, then we *must* split the calculation
0510    // as the size of the ULP changes with each order of magnitude change:
0511    //
0512    if(b > upper)
0513    {
0514       int expon2;
0515       (void)frexp(b, &expon2);
0516       T upper2 = ldexp(T(0.5), expon2);
0517       result = float_distance(upper2, b);
0518       result += (expon2 - expon - 1) * ldexp(T(1), tools::digits<T>() - 1);
0519    }
0520    //
0521    // Use compensated double-double addition to avoid rounding
0522    // errors in the subtraction:
0523    //
0524    expon = tools::digits<T>() - expon;
0525    T mb, x, y, z;
0526    if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
0527    {
0528       //
0529       // Special case - either one end of the range is a denormal, or else the difference is.
0530       // The regular code will fail if we're using the SSE2 registers on Intel and either
0531       // the FTZ or DAZ flags are set.
0532       //
0533       T a2 = ldexp(a, tools::digits<T>());
0534       T b2 = ldexp(b, tools::digits<T>());
0535       mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2);
0536       x = a2 + mb;
0537       z = x - a2;
0538       y = (a2 - (x - z)) + (mb - z);
0539
0540       expon -= tools::digits<T>();
0541    }
0542    else
0543    {
0544       mb = -(std::min)(upper, b);
0545       x = a + mb;
0546       z = x - a;
0547       y = (a - (x - z)) + (mb - z);
0548    }
0549    if(x < 0)
0550    {
0551       x = -x;
0552       y = -y;
0553    }
0554    result += ldexp(x, expon) + ldexp(y, expon);
0555    //
0556    // Result must be an integer:
0557    //
0558    BOOST_MATH_ASSERT(result == floor(result));
0559    return result;
0560 } // float_distance_imp
0561 //
0562 // Special versions for bases other than 2:
0563 //
0564 template <class T, class Policy>
0565 T float_distance_imp(const T& a, const T& b, const std::false_type&, const Policy& pol)
0566 {
0567    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0568    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0569
0570    BOOST_MATH_STD_USING
0571    //
0572    // Error handling:
0573    //
0574    static const char* function = "float_distance<%1%>(%1%, %1%)";
0575    if(!(boost::math::isfinite)(a))
0576       return policies::raise_domain_error<T>(function, "Argument a must be finite, but got %1%", a, pol);
0577    if(!(boost::math::isfinite)(b))
0578       return policies::raise_domain_error<T>(function, "Argument b must be finite, but got %1%", b, pol);
0579    //
0580    // Special cases:
0581    //
0582    if(a > b)
0583       return -float_distance(b, a, pol);
0584    if(a == b)
0585       return T(0);
0586    if(a == 0)
0587       return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
0588    if(b == 0)
0589       return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
0590    if(boost::math::sign(a) != boost::math::sign(b))
0591       return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
0592          + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
0593    //
0594    // By the time we get here, both a and b must have the same sign, we want
0595    // b > a and both positive for the following logic:
0596    //
0597    if(a < 0)
0598       return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
0599
0600    BOOST_MATH_ASSERT(a >= 0);
0601    BOOST_MATH_ASSERT(b >= a);
0602
0603    std::intmax_t expon;
0604    //
0605    // Note that if a is a denorm then the usual formula fails
0606    // because we actually have fewer than tools::digits<T>()
0607    // significant bits in the representation:
0608    //
0609    expon = 1 + ilogb(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a);
0610    T upper = scalbn(T(1), expon);
0611    T result = T(0);
0612    //
0613    // If b is greater than upper, then we *must* split the calculation
0614    // as the size of the ULP changes with each order of magnitude change:
0615    //
0616    if(b > upper)
0617    {
0618       std::intmax_t expon2 = 1 + ilogb(b);
0619       T upper2 = scalbn(T(1), expon2 - 1);
0620       result = float_distance(upper2, b);
0621       result += (expon2 - expon - 1) * scalbn(T(1), std::numeric_limits<T>::digits - 1);
0622    }
0623    //
0624    // Use compensated double-double addition to avoid rounding
0625    // errors in the subtraction:
0626    //
0627    expon = std::numeric_limits<T>::digits - expon;
0628    T mb, x, y, z;
0629    if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
0630    {
0631       //
0632       // Special case - either one end of the range is a denormal, or else the difference is.
0633       // The regular code will fail if we're using the SSE2 registers on Intel and either
0634       // the FTZ or DAZ flags are set.
0635       //
0636       T a2 = scalbn(a, std::numeric_limits<T>::digits);
0637       T b2 = scalbn(b, std::numeric_limits<T>::digits);
0638       mb = -(std::min)(T(scalbn(upper, std::numeric_limits<T>::digits)), b2);
0639       x = a2 + mb;
0640       z = x - a2;
0641       y = (a2 - (x - z)) + (mb - z);
0642
0643       expon -= std::numeric_limits<T>::digits;
0644    }
0645    else
0646    {
0647       mb = -(std::min)(upper, b);
0648       x = a + mb;
0649       z = x - a;
0650       y = (a - (x - z)) + (mb - z);
0651    }
0652    if(x < 0)
0653    {
0654       x = -x;
0655       y = -y;
0656    }
0657    result += scalbn(x, expon) + scalbn(y, expon);
0658    //
0659    // Result must be an integer:
0660    //
0661    BOOST_MATH_ASSERT(result == floor(result));
0662    return result;
0663 } // float_distance_imp
0664
0665 } // namespace detail
0666
0667 template <class T, class U, class Policy>
0668 inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
0669 {
0670    //
0671    // We allow ONE of a and b to be an integer type, otherwise both must be the SAME type.
0672    //
0673    static_assert(
0674       (std::is_same<T, U>::value
0675       || (std::is_integral<T>::value && !std::is_integral<U>::value)
0676       || (!std::is_integral<T>::value && std::is_integral<U>::value)
0677       || (std::numeric_limits<T>::is_specialized && std::numeric_limits<U>::is_specialized
0678          && (std::numeric_limits<T>::digits == std::numeric_limits<U>::digits)
0679          && (std::numeric_limits<T>::radix == std::numeric_limits<U>::radix)
0680          && !std::numeric_limits<T>::is_integer && !std::numeric_limits<U>::is_integer)),
0681       "Float distance between two different floating point types is undefined.");
0682
0683    BOOST_MATH_IF_CONSTEXPR (!std::is_same<T, U>::value)
0684    {
0685       BOOST_MATH_IF_CONSTEXPR(std::is_integral<T>::value)
0686       {
0687          return float_distance(static_cast<U>(a), b, pol);
0688       }
0689       else
0690       {
0691          return float_distance(a, static_cast<T>(b), pol);
0692       }
0693    }
0694    else
0695    {
0696       typedef typename tools::promote_args<T, U>::type result_type;
0697       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()), std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
0698    }
0699 }
0700
0701 template <class T, class U>
0702 typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
0703 {
0704    return boost::math::float_distance(a, b, policies::policy<>());
0705 }
0706
0707 namespace detail{
0708
0709 template <class T, class Policy>
0710 T float_advance_imp(T val, int distance, const std::true_type&, const Policy& pol)
0711 {
0712    BOOST_MATH_STD_USING
0713    //
0714    // Error handling:
0715    //
0716    static const char* function = "float_advance<%1%>(%1%, int)";
0717
0718    int fpclass = (boost::math::fpclassify)(val);
0719
0720    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
0721       return policies::raise_domain_error<T>(function, "Argument val must be finite, but got %1%", val, pol);
0722
0723    if(val < 0)
0724       return -float_advance(-val, -distance, pol);
0725    if(distance == 0)
0726       return val;
0727    if(distance == 1)
0728       return float_next(val, pol);
0729    if(distance == -1)
0730       return float_prior(val, pol);
0731
0732    if(fabs(val) < detail::get_min_shift_value<T>())
0733    {
0734       //
0735       // Special case: if the value of the least significant bit is a denorm,
0736       // implement in terms of float_next/float_prior.
0737       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0738       //
0739       if(distance > 0)
0740       {
0741          do{ val = float_next(val, pol); } while(--distance);
0742       }
0743       else
0744       {
0745          do{ val = float_prior(val, pol); } while(++distance);
0746       }
0747       return val;
0748    }
0749
0750    int expon;
0751    (void)frexp(val, &expon);
0752    T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon);
0753    // We can not have denorms here, since we have taken care of them above:
0754    BOOST_MATH_ASSERT(val > tools::min_value<T>());
0755    T limit_distance = float_distance(val, limit);
0756    while(fabs(limit_distance) < abs(distance))
0757    {
0758       distance -= itrunc(limit_distance);
0759       val = limit;
0760       if(distance < 0)
0761       {
0762          limit /= 2;
0763          expon--;
0764       }
0765       else
0766       {
0767          limit *= 2;
0768          expon++;
0769       }
0770       limit_distance = float_distance(val, limit);
0771       if(distance && (limit_distance == 0))
0772       {
0773          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);  // LCOV_EXCL_LINE This *should* be unreachable.
0774       }
0775    }
0776    if((0.5f == frexp(val, &expon)) && (distance < 0))
0777       --expon;
0778    T diff = 0;
0779    if(val != 0)
0780       diff = distance * ldexp(T(1), expon - tools::digits<T>());
0781    if(diff == 0)
0782       diff = distance * detail::get_smallest_value<T>(); // LCOV_EXCL_LINE This *should* be unreachable given that denorms are handled above already.
0783    return val += diff;
0784 } // float_advance_imp
0785 //
0786 // Special version for bases other than 2:
0787 //
0788 template <class T, class Policy>
0789 T float_advance_imp(T val, int distance, const std::false_type&, const Policy& pol)
0790 {
0791    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0792    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0793
0794    BOOST_MATH_STD_USING
0795    //
0796    // Error handling:
0797    //
0798    static const char* function = "float_advance<%1%>(%1%, int)";
0799
0800    int fpclass = (boost::math::fpclassify)(val);
0801
0802    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
0803       return policies::raise_domain_error<T>(function, "Argument val must be finite, but got %1%", val, pol);
0804
0805    if(val < 0)
0806       return -float_advance(-val, -distance, pol);
0807    if(distance == 0)
0808       return val;
0809    if(distance == 1)
0810       return float_next(val, pol);
0811    if(distance == -1)
0812       return float_prior(val, pol);
0813
0814    if(fabs(val) < detail::get_min_shift_value<T>())
0815    {
0816       //
0817       // Special case: if the value of the least significant bit is a denorm,
0818       // implement in terms of float_next/float_prior.
0819       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0820       //
0821       if(distance > 0)
0822       {
0823          do{ val = float_next(val, pol); } while(--distance);
0824       }
0825       else
0826       {
0827          do{ val = float_prior(val, pol); } while(++distance);
0828       }
0829       return val;
0830    }
0831
0832    std::intmax_t expon = 1 + ilogb(val);
0833    T limit = scalbn(T(1), distance < 0 ? expon - 1 : expon);
0834    BOOST_MATH_ASSERT(val > tools::min_value<T>()); // denorms already handled.
0835    T limit_distance = float_distance(val, limit);
0836    while(fabs(limit_distance) < abs(distance))
0837    {
0838       distance -= itrunc(limit_distance);
0839       val = limit;
0840       if(distance < 0)
0841       {
0842          limit /= std::numeric_limits<T>::radix;
0843          expon--;
0844       }
0845       else
0846       {
0847          limit *= std::numeric_limits<T>::radix; // LCOV_EXCL_LINE Probably unreachable for the decimal types we have?
0848          expon++;                                // LCOV_EXCL_LINE
0849       }
0850       limit_distance = float_distance(val, limit);
0851       if(distance && (limit_distance == 0))
0852       {
0853          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);  // LCOV_EXCL_LINE should never get here!
0854       }
0855    }
0856    /*expon = 1 + ilogb(val);
0857    if((1 == scalbn(val, 1 + expon)) && (distance < 0))
0858       --expon;*/
0859    T diff = 0;
0860    if(val != 0)
0861       diff = distance * scalbn(T(1), expon - std::numeric_limits<T>::digits);
0862    if(diff == 0)
0863       diff = distance * detail::get_smallest_value<T>(); // LCOV_EXCL_LINE This *should* be unreachable given that denorms are handled above.
0864    return val += diff;
0865 } // float_advance_imp
0866
0867 } // namespace detail
0868
0869 template <class T, class Policy>
0870 inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
0871 {
0872    typedef typename tools::promote_args<T>::type result_type;
0873    return detail::float_advance_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), distance, std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
0874 }
0875
0876 template <class T>
0877 inline typename tools::promote_args<T>::type float_advance(const T& val, int distance)
0878 {
0879    return boost::math::float_advance(val, distance, policies::policy<>());
0880 }
0881
0882 }} // boost math namespaces
0883
0884 #endif
0885
0886 #endif // BOOST_MATH_SPECIAL_NEXT_HPP