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