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