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 struct min_shift_initializer
0131 {
0132    struct init
0133    {
0134       init()
0135       {
0136          do_init();
0137       }
0138       static void do_init()
0139       {
0140          get_min_shift_value<T>();
0141       }
0142       void force_instantiate()const{}
0143    };
0144    static const init initializer;
0145    static void force_instantiate()
0146    {
0147       initializer.force_instantiate();
0148    }
0149 };
0150
0151 template <class T>
0152 const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
0153
0154 template <class T>
0155 inline T calc_min_shifted(const std::true_type&)
0156 {
0157    BOOST_MATH_STD_USING
0158    return ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
0159 }
0160 template <class T>
0161 inline T calc_min_shifted(const std::false_type&)
0162 {
0163    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0164    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0165
0166    return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);
0167 }
0168
0169
0170 template <class T>
0171 inline T get_min_shift_value()
0172 {
0173    static const T val = calc_min_shifted<T>(std::integral_constant<bool, !std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());
0174    min_shift_initializer<T>::force_instantiate();
0175
0176    return val;
0177 }
0178
0179 template <class T, bool b = boost::math::tools::detail::has_backend_type<T>::value>
0180 struct exponent_type
0181 {
0182    typedef int type;
0183 };
0184
0185 template <class T>
0186 struct exponent_type<T, true>
0187 {
0188    typedef typename T::backend_type::exponent_type type;
0189 };
0190
0191 template <class T, class Policy>
0192 T float_next_imp(const T& val, const std::true_type&, const Policy& pol)
0193 {
0194    typedef typename exponent_type<T>::type exponent_type;
0195
0196    BOOST_MATH_STD_USING
0197    exponent_type expon;
0198    static const char* function = "float_next<%1%>(%1%)";
0199
0200    int fpclass = (boost::math::fpclassify)(val);
0201
0202    if (fpclass == (int)FP_INFINITE)
0203    {
0204       if (val < 0)
0205          return -tools::max_value<T>();
0206       return val;  // +INF
0207    }
0208    else if (fpclass == (int)FP_NAN)
0209    {
0210       return policies::raise_domain_error<T>(
0211          function,
0212          "Argument must be finite, but got %1%", val, pol);
0213    }
0214
0215    if(val >= tools::max_value<T>())
0216       return policies::raise_overflow_error<T>(function, nullptr, pol);
0217
0218    if(val == 0)
0219       return detail::get_smallest_value<T>();
0220
0221    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
0222    {
0223       //
0224       // Special case: if the value of the least significant bit is a denorm, and the result
0225       // would not be a denorm, then shift the input, increment, and shift back.
0226       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0227       //
0228       return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
0229    }
0230
0231    if(-0.5f == frexp(val, &expon))
0232       --expon; // reduce exponent when val is a power of two, and negative.
0233    T diff = ldexp(T(1), expon - tools::digits<T>());
0234    if(diff == 0)
0235       diff = detail::get_smallest_value<T>();
0236    return val + diff;
0237 } // float_next_imp
0238 //
0239 // Special version for some base other than 2:
0240 //
0241 template <class T, class Policy>
0242 T float_next_imp(const T& val, const std::false_type&, const Policy& pol)
0243 {
0244    typedef typename exponent_type<T>::type exponent_type;
0245
0246    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0247    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0248
0249    BOOST_MATH_STD_USING
0250    exponent_type expon;
0251    static const char* function = "float_next<%1%>(%1%)";
0252
0253    int fpclass = (boost::math::fpclassify)(val);
0254
0255    if (fpclass == (int)FP_INFINITE)
0256    {
0257       if (val < 0)
0258          return -tools::max_value<T>();
0259       return val;  // +INF
0260    }
0261    else if (fpclass == (int)FP_NAN)
0262    {
0263       return policies::raise_domain_error<T>(
0264          function,
0265          "Argument must be finite, but got %1%", val, pol);
0266    }
0267
0268    if(val >= tools::max_value<T>())
0269       return policies::raise_overflow_error<T>(function, nullptr, pol);
0270
0271    if(val == 0)
0272       return detail::get_smallest_value<T>();
0273
0274    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
0275    {
0276       //
0277       // Special case: if the value of the least significant bit is a denorm, and the result
0278       // would not be a denorm, then shift the input, increment, and shift back.
0279       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0280       //
0281       return scalbn(float_next(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
0282    }
0283
0284    expon = 1 + ilogb(val);
0285    if(-1 == scalbn(val, -expon) * std::numeric_limits<T>::radix)
0286       --expon; // reduce exponent when val is a power of base, and negative.
0287    T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
0288    if(diff == 0)
0289       diff = detail::get_smallest_value<T>();
0290    return val + diff;
0291 } // float_next_imp
0292
0293 } // namespace detail
0294
0295 template <class T, class Policy>
0296 inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
0297 {
0298    typedef typename tools::promote_args<T>::type result_type;
0299    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);
0300 }
0301
0302 #if 0 //def BOOST_MSVC
0303 //
0304 // We used to use ::_nextafter here, but doing so fails when using
0305 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
0306 // - albeit slower - code instead as at least that gives the correct answer.
0307 //
0308 template <class Policy>
0309 inline double float_next(const double& val, const Policy& pol)
0310 {
0311    static const char* function = "float_next<%1%>(%1%)";
0312
0313    if(!(boost::math::isfinite)(val) && (val > 0))
0314       return policies::raise_domain_error<double>(
0315          function,
0316          "Argument must be finite, but got %1%", val, pol);
0317
0318    if(val >= tools::max_value<double>())
0319       return policies::raise_overflow_error<double>(function, nullptr, pol);
0320
0321    return ::_nextafter(val, tools::max_value<double>());
0322 }
0323 #endif
0324
0325 template <class T>
0326 inline typename tools::promote_args<T>::type float_next(const T& val)
0327 {
0328    return float_next(val, policies::policy<>());
0329 }
0330
0331 namespace detail{
0332
0333 template <class T, class Policy>
0334 T float_prior_imp(const T& val, const std::true_type&, const Policy& pol)
0335 {
0336    typedef typename exponent_type<T>::type exponent_type;
0337
0338    BOOST_MATH_STD_USING
0339    exponent_type expon;
0340    static const char* function = "float_prior<%1%>(%1%)";
0341
0342    int fpclass = (boost::math::fpclassify)(val);
0343
0344    if (fpclass == (int)FP_INFINITE)
0345    {
0346       if (val > 0)
0347          return tools::max_value<T>();
0348       return val; // -INF
0349    }
0350    else if (fpclass == (int)FP_NAN)
0351    {
0352       return policies::raise_domain_error<T>(
0353          function,
0354          "Argument must be finite, but got %1%", val, pol);
0355    }
0356
0357    if(val <= -tools::max_value<T>())
0358       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0359
0360    if(val == 0)
0361       return -detail::get_smallest_value<T>();
0362
0363    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
0364    {
0365       //
0366       // Special case: if the value of the least significant bit is a denorm, and the result
0367       // would not be a denorm, then shift the input, increment, and shift back.
0368       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0369       //
0370       return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
0371    }
0372
0373    T remain = frexp(val, &expon);
0374    if(remain == 0.5f)
0375       --expon; // when val is a power of two we must reduce the exponent
0376    T diff = ldexp(T(1), expon - tools::digits<T>());
0377    if(diff == 0)
0378       diff = detail::get_smallest_value<T>();
0379    return val - diff;
0380 } // float_prior_imp
0381 //
0382 // Special version for bases other than 2:
0383 //
0384 template <class T, class Policy>
0385 T float_prior_imp(const T& val, const std::false_type&, const Policy& pol)
0386 {
0387    typedef typename exponent_type<T>::type exponent_type;
0388
0389    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0390    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0391
0392    BOOST_MATH_STD_USING
0393    exponent_type expon;
0394    static const char* function = "float_prior<%1%>(%1%)";
0395
0396    int fpclass = (boost::math::fpclassify)(val);
0397
0398    if (fpclass == (int)FP_INFINITE)
0399    {
0400       if (val > 0)
0401          return tools::max_value<T>();
0402       return val; // -INF
0403    }
0404    else if (fpclass == (int)FP_NAN)
0405    {
0406       return policies::raise_domain_error<T>(
0407          function,
0408          "Argument must be finite, but got %1%", val, pol);
0409    }
0410
0411    if(val <= -tools::max_value<T>())
0412       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0413
0414    if(val == 0)
0415       return -detail::get_smallest_value<T>();
0416
0417    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
0418    {
0419       //
0420       // Special case: if the value of the least significant bit is a denorm, and the result
0421       // would not be a denorm, then shift the input, increment, and shift back.
0422       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0423       //
0424       return scalbn(float_prior(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
0425    }
0426
0427    expon = 1 + ilogb(val);
0428    T remain = scalbn(val, -expon);
0429    if(remain * std::numeric_limits<T>::radix == 1)
0430       --expon; // when val is a power of two we must reduce the exponent
0431    T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
0432    if(diff == 0)
0433       diff = detail::get_smallest_value<T>();
0434    return val - diff;
0435 } // float_prior_imp
0436
0437 } // namespace detail
0438
0439 template <class T, class Policy>
0440 inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
0441 {
0442    typedef typename tools::promote_args<T>::type result_type;
0443    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);
0444 }
0445
0446 #if 0 //def BOOST_MSVC
0447 //
0448 // We used to use ::_nextafter here, but doing so fails when using
0449 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
0450 // - albeit slower - code instead as at least that gives the correct answer.
0451 //
0452 template <class Policy>
0453 inline double float_prior(const double& val, const Policy& pol)
0454 {
0455    static const char* function = "float_prior<%1%>(%1%)";
0456
0457    if(!(boost::math::isfinite)(val) && (val < 0))
0458       return policies::raise_domain_error<double>(
0459          function,
0460          "Argument must be finite, but got %1%", val, pol);
0461
0462    if(val <= -tools::max_value<double>())
0463       return -policies::raise_overflow_error<double>(function, nullptr, pol);
0464
0465    return ::_nextafter(val, -tools::max_value<double>());
0466 }
0467 #endif
0468
0469 template <class T>
0470 inline typename tools::promote_args<T>::type float_prior(const T& val)
0471 {
0472    return float_prior(val, policies::policy<>());
0473 }
0474
0475 template <class T, class U, class Policy>
0476 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol)
0477 {
0478    typedef typename tools::promote_args<T, U>::type result_type;
0479    return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol);
0480 }
0481
0482 template <class T, class U>
0483 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
0484 {
0485    return nextafter(val, direction, policies::policy<>());
0486 }
0487
0488 namespace detail{
0489
0490 template <class T, class Policy>
0491 T float_distance_imp(const T& a, const T& b, const std::true_type&, const Policy& pol)
0492 {
0493    BOOST_MATH_STD_USING
0494    //
0495    // Error handling:
0496    //
0497    static const char* function = "float_distance<%1%>(%1%, %1%)";
0498    if(!(boost::math::isfinite)(a))
0499       return policies::raise_domain_error<T>(
0500          function,
0501          "Argument a must be finite, but got %1%", a, pol);
0502    if(!(boost::math::isfinite)(b))
0503       return policies::raise_domain_error<T>(
0504          function,
0505          "Argument b must be finite, but got %1%", b, pol);
0506    //
0507    // Special cases:
0508    //
0509    if(a > b)
0510       return -float_distance(b, a, pol);
0511    if(a == b)
0512       return T(0);
0513    if(a == 0)
0514       return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
0515    if(b == 0)
0516       return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
0517    if(boost::math::sign(a) != boost::math::sign(b))
0518       return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
0519          + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
0520    //
0521    // By the time we get here, both a and b must have the same sign, we want
0522    // b > a and both positive for the following logic:
0523    //
0524    if(a < 0)
0525       return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
0526
0527    BOOST_MATH_ASSERT(a >= 0);
0528    BOOST_MATH_ASSERT(b >= a);
0529
0530    int expon;
0531    //
0532    // Note that if a is a denorm then the usual formula fails
0533    // because we actually have fewer than tools::digits<T>()
0534    // significant bits in the representation:
0535    //
0536    (void)frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
0537    T upper = ldexp(T(1), expon);
0538    T result = T(0);
0539    //
0540    // If b is greater than upper, then we *must* split the calculation
0541    // as the size of the ULP changes with each order of magnitude change:
0542    //
0543    if(b > upper)
0544    {
0545       int expon2;
0546       (void)frexp(b, &expon2);
0547       T upper2 = ldexp(T(0.5), expon2);
0548       result = float_distance(upper2, b);
0549       result += (expon2 - expon - 1) * ldexp(T(1), tools::digits<T>() - 1);
0550    }
0551    //
0552    // Use compensated double-double addition to avoid rounding
0553    // errors in the subtraction:
0554    //
0555    expon = tools::digits<T>() - expon;
0556    T mb, x, y, z;
0557    if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
0558    {
0559       //
0560       // Special case - either one end of the range is a denormal, or else the difference is.
0561       // The regular code will fail if we're using the SSE2 registers on Intel and either
0562       // the FTZ or DAZ flags are set.
0563       //
0564       T a2 = ldexp(a, tools::digits<T>());
0565       T b2 = ldexp(b, tools::digits<T>());
0566       mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2);
0567       x = a2 + mb;
0568       z = x - a2;
0569       y = (a2 - (x - z)) + (mb - z);
0570
0571       expon -= tools::digits<T>();
0572    }
0573    else
0574    {
0575       mb = -(std::min)(upper, b);
0576       x = a + mb;
0577       z = x - a;
0578       y = (a - (x - z)) + (mb - z);
0579    }
0580    if(x < 0)
0581    {
0582       x = -x;
0583       y = -y;
0584    }
0585    result += ldexp(x, expon) + ldexp(y, expon);
0586    //
0587    // Result must be an integer:
0588    //
0589    BOOST_MATH_ASSERT(result == floor(result));
0590    return result;
0591 } // float_distance_imp
0592 //
0593 // Special versions for bases other than 2:
0594 //
0595 template <class T, class Policy>
0596 T float_distance_imp(const T& a, const T& b, const std::false_type&, const Policy& pol)
0597 {
0598    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0599    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0600
0601    BOOST_MATH_STD_USING
0602    //
0603    // Error handling:
0604    //
0605    static const char* function = "float_distance<%1%>(%1%, %1%)";
0606    if(!(boost::math::isfinite)(a))
0607       return policies::raise_domain_error<T>(
0608          function,
0609          "Argument a must be finite, but got %1%", a, pol);
0610    if(!(boost::math::isfinite)(b))
0611       return policies::raise_domain_error<T>(
0612          function,
0613          "Argument b must be finite, but got %1%", b, pol);
0614    //
0615    // Special cases:
0616    //
0617    if(a > b)
0618       return -float_distance(b, a, pol);
0619    if(a == b)
0620       return T(0);
0621    if(a == 0)
0622       return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
0623    if(b == 0)
0624       return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
0625    if(boost::math::sign(a) != boost::math::sign(b))
0626       return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
0627          + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
0628    //
0629    // By the time we get here, both a and b must have the same sign, we want
0630    // b > a and both positive for the following logic:
0631    //
0632    if(a < 0)
0633       return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
0634
0635    BOOST_MATH_ASSERT(a >= 0);
0636    BOOST_MATH_ASSERT(b >= a);
0637
0638    std::intmax_t expon;
0639    //
0640    // Note that if a is a denorm then the usual formula fails
0641    // because we actually have fewer than tools::digits<T>()
0642    // significant bits in the representation:
0643    //
0644    expon = 1 + ilogb(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a);
0645    T upper = scalbn(T(1), expon);
0646    T result = T(0);
0647    //
0648    // If b is greater than upper, then we *must* split the calculation
0649    // as the size of the ULP changes with each order of magnitude change:
0650    //
0651    if(b > upper)
0652    {
0653       std::intmax_t expon2 = 1 + ilogb(b);
0654       T upper2 = scalbn(T(1), expon2 - 1);
0655       result = float_distance(upper2, b);
0656       result += (expon2 - expon - 1) * scalbn(T(1), std::numeric_limits<T>::digits - 1);
0657    }
0658    //
0659    // Use compensated double-double addition to avoid rounding
0660    // errors in the subtraction:
0661    //
0662    expon = std::numeric_limits<T>::digits - expon;
0663    T mb, x, y, z;
0664    if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
0665    {
0666       //
0667       // Special case - either one end of the range is a denormal, or else the difference is.
0668       // The regular code will fail if we're using the SSE2 registers on Intel and either
0669       // the FTZ or DAZ flags are set.
0670       //
0671       T a2 = scalbn(a, std::numeric_limits<T>::digits);
0672       T b2 = scalbn(b, std::numeric_limits<T>::digits);
0673       mb = -(std::min)(T(scalbn(upper, std::numeric_limits<T>::digits)), b2);
0674       x = a2 + mb;
0675       z = x - a2;
0676       y = (a2 - (x - z)) + (mb - z);
0677
0678       expon -= std::numeric_limits<T>::digits;
0679    }
0680    else
0681    {
0682       mb = -(std::min)(upper, b);
0683       x = a + mb;
0684       z = x - a;
0685       y = (a - (x - z)) + (mb - z);
0686    }
0687    if(x < 0)
0688    {
0689       x = -x;
0690       y = -y;
0691    }
0692    result += scalbn(x, expon) + scalbn(y, expon);
0693    //
0694    // Result must be an integer:
0695    //
0696    BOOST_MATH_ASSERT(result == floor(result));
0697    return result;
0698 } // float_distance_imp
0699
0700 } // namespace detail
0701
0702 template <class T, class U, class Policy>
0703 inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
0704 {
0705    //
0706    // We allow ONE of a and b to be an integer type, otherwise both must be the SAME type.
0707    //
0708    static_assert(
0709       (std::is_same<T, U>::value
0710       || (std::is_integral<T>::value && !std::is_integral<U>::value)
0711       || (!std::is_integral<T>::value && std::is_integral<U>::value)
0712       || (std::numeric_limits<T>::is_specialized && std::numeric_limits<U>::is_specialized
0713          && (std::numeric_limits<T>::digits == std::numeric_limits<U>::digits)
0714          && (std::numeric_limits<T>::radix == std::numeric_limits<U>::radix)
0715          && !std::numeric_limits<T>::is_integer && !std::numeric_limits<U>::is_integer)),
0716       "Float distance between two different floating point types is undefined.");
0717
0718    BOOST_MATH_IF_CONSTEXPR (!std::is_same<T, U>::value)
0719    {
0720       BOOST_MATH_IF_CONSTEXPR(std::is_integral<T>::value)
0721       {
0722          return float_distance(static_cast<U>(a), b, pol);
0723       }
0724       else
0725       {
0726          return float_distance(a, static_cast<T>(b), pol);
0727       }
0728    }
0729    else
0730    {
0731       typedef typename tools::promote_args<T, U>::type result_type;
0732       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);
0733    }
0734 }
0735
0736 template <class T, class U>
0737 typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
0738 {
0739    return boost::math::float_distance(a, b, policies::policy<>());
0740 }
0741
0742 namespace detail{
0743
0744 template <class T, class Policy>
0745 T float_advance_imp(T val, int distance, const std::true_type&, const Policy& pol)
0746 {
0747    BOOST_MATH_STD_USING
0748    //
0749    // Error handling:
0750    //
0751    static const char* function = "float_advance<%1%>(%1%, int)";
0752
0753    int fpclass = (boost::math::fpclassify)(val);
0754
0755    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
0756       return policies::raise_domain_error<T>(
0757          function,
0758          "Argument val must be finite, but got %1%", val, pol);
0759
0760    if(val < 0)
0761       return -float_advance(-val, -distance, pol);
0762    if(distance == 0)
0763       return val;
0764    if(distance == 1)
0765       return float_next(val, pol);
0766    if(distance == -1)
0767       return float_prior(val, pol);
0768
0769    if(fabs(val) < detail::get_min_shift_value<T>())
0770    {
0771       //
0772       // Special case: if the value of the least significant bit is a denorm,
0773       // implement in terms of float_next/float_prior.
0774       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0775       //
0776       if(distance > 0)
0777       {
0778          do{ val = float_next(val, pol); } while(--distance);
0779       }
0780       else
0781       {
0782          do{ val = float_prior(val, pol); } while(++distance);
0783       }
0784       return val;
0785    }
0786
0787    int expon;
0788    (void)frexp(val, &expon);
0789    T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon);
0790    if(val <= tools::min_value<T>())
0791    {
0792       limit = sign(T(distance)) * tools::min_value<T>();
0793    }
0794    T limit_distance = float_distance(val, limit);
0795    while(fabs(limit_distance) < abs(distance))
0796    {
0797       distance -= itrunc(limit_distance);
0798       val = limit;
0799       if(distance < 0)
0800       {
0801          limit /= 2;
0802          expon--;
0803       }
0804       else
0805       {
0806          limit *= 2;
0807          expon++;
0808       }
0809       limit_distance = float_distance(val, limit);
0810       if(distance && (limit_distance == 0))
0811       {
0812          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);
0813       }
0814    }
0815    if((0.5f == frexp(val, &expon)) && (distance < 0))
0816       --expon;
0817    T diff = 0;
0818    if(val != 0)
0819       diff = distance * ldexp(T(1), expon - tools::digits<T>());
0820    if(diff == 0)
0821       diff = distance * detail::get_smallest_value<T>();
0822    return val += diff;
0823 } // float_advance_imp
0824 //
0825 // Special version for bases other than 2:
0826 //
0827 template <class T, class Policy>
0828 T float_advance_imp(T val, int distance, const std::false_type&, const Policy& pol)
0829 {
0830    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0831    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0832
0833    BOOST_MATH_STD_USING
0834    //
0835    // Error handling:
0836    //
0837    static const char* function = "float_advance<%1%>(%1%, int)";
0838
0839    int fpclass = (boost::math::fpclassify)(val);
0840
0841    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
0842       return policies::raise_domain_error<T>(
0843          function,
0844          "Argument val must be finite, but got %1%", val, pol);
0845
0846    if(val < 0)
0847       return -float_advance(-val, -distance, pol);
0848    if(distance == 0)
0849       return val;
0850    if(distance == 1)
0851       return float_next(val, pol);
0852    if(distance == -1)
0853       return float_prior(val, pol);
0854
0855    if(fabs(val) < detail::get_min_shift_value<T>())
0856    {
0857       //
0858       // Special case: if the value of the least significant bit is a denorm,
0859       // implement in terms of float_next/float_prior.
0860       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0861       //
0862       if(distance > 0)
0863       {
0864          do{ val = float_next(val, pol); } while(--distance);
0865       }
0866       else
0867       {
0868          do{ val = float_prior(val, pol); } while(++distance);
0869       }
0870       return val;
0871    }
0872
0873    std::intmax_t expon = 1 + ilogb(val);
0874    T limit = scalbn(T(1), distance < 0 ? expon - 1 : expon);
0875    if(val <= tools::min_value<T>())
0876    {
0877       limit = sign(T(distance)) * tools::min_value<T>();
0878    }
0879    T limit_distance = float_distance(val, limit);
0880    while(fabs(limit_distance) < abs(distance))
0881    {
0882       distance -= itrunc(limit_distance);
0883       val = limit;
0884       if(distance < 0)
0885       {
0886          limit /= std::numeric_limits<T>::radix;
0887          expon--;
0888       }
0889       else
0890       {
0891          limit *= std::numeric_limits<T>::radix;
0892          expon++;
0893       }
0894       limit_distance = float_distance(val, limit);
0895       if(distance && (limit_distance == 0))
0896       {
0897          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);
0898       }
0899    }
0900    /*expon = 1 + ilogb(val);
0901    if((1 == scalbn(val, 1 + expon)) && (distance < 0))
0902       --expon;*/
0903    T diff = 0;
0904    if(val != 0)
0905       diff = distance * scalbn(T(1), expon - std::numeric_limits<T>::digits);
0906    if(diff == 0)
0907       diff = distance * detail::get_smallest_value<T>();
0908    return val += diff;
0909 } // float_advance_imp
0910
0911 } // namespace detail
0912
0913 template <class T, class Policy>
0914 inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
0915 {
0916    typedef typename tools::promote_args<T>::type result_type;
0917    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);
0918 }
0919
0920 template <class T>
0921 inline typename tools::promote_args<T>::type float_advance(const T& val, int distance)
0922 {
0923    return boost::math::float_advance(val, distance, policies::policy<>());
0924 }
0925
0926 }} // boost math namespaces
0927
0928 #endif
0929
0930 #endif // BOOST_MATH_SPECIAL_NEXT_HPP