boost/math/special_functions/next.hpp

   1 //  (C) Copyright John Maddock 2008.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_SPECIAL_NEXT_HPP
   7 #define BOOST_MATH_SPECIAL_NEXT_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/special_functions/math_fwd.hpp>
  14 #include <boost/math/policies/error_handling.hpp>
  15 #include <boost/math/special_functions/fpclassify.hpp>
  16 #include <boost/math/special_functions/sign.hpp>
  17 #include <boost/math/special_functions/trunc.hpp>
  18
  19 #include <float.h>
  20
  21 #if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3)))
  22 #if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__)
  23 #include "xmmintrin.h"
  24 #define BOOST_MATH_CHECK_SSE2
  25 #endif
  26 #endif
  27
  28 namespace boost{ namespace math{
  29
  30    namespace concepts {
  31
  32       class real_concept;
  33       class std_real_concept;
  34
  35    }
  36
  37 namespace detail{
  38
  39 template <class T>
  40 struct has_hidden_guard_digits;
  41 template <>
  42 struct has_hidden_guard_digits<float> : public mpl::false_ {};
  43 template <>
  44 struct has_hidden_guard_digits<double> : public mpl::false_ {};
  45 template <>
  46 struct has_hidden_guard_digits<long double> : public mpl::false_ {};
  47 #ifdef BOOST_HAS_FLOAT128
  48 template <>
  49 struct has_hidden_guard_digits<__float128> : public mpl::false_ {};
  50 #endif
  51 template <>
  52 struct has_hidden_guard_digits<boost::math::concepts::real_concept> : public mpl::false_ {};
  53 template <>
  54 struct has_hidden_guard_digits<boost::math::concepts::std_real_concept> : public mpl::false_ {};
  55
  56 template <class T, bool b>
  57 struct has_hidden_guard_digits_10 : public mpl::false_ {};
  58 template <class T>
  59 struct has_hidden_guard_digits_10<T, true> : public mpl::bool_<(std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {};
  60
  61 template <class T>
  62 struct has_hidden_guard_digits
  63    : public has_hidden_guard_digits_10<T,
  64    std::numeric_limits<T>::is_specialized
  65    && (std::numeric_limits<T>::radix == 10) >
  66 {};
  67
  68 template <class T>
  69 inline const T& normalize_value(const T& val, const mpl::false_&) { return val; }
  70 template <class T>
  71 inline T normalize_value(const T& val, const mpl::true_&)
  72 {
  73    BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
  74    BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
  75
  76    boost::intmax_t shift = (boost::intmax_t)std::numeric_limits<T>::digits - (boost::intmax_t)ilogb(val) - 1;
  77    T result = scalbn(val, shift);
  78    result = round(result);
  79    return scalbn(result, -shift);
  80 }
  81
  82 template <class T>
  83 inline T get_smallest_value(mpl::true_ const&)
  84 {
  85    //
  86    // numeric_limits lies about denorms being present - particularly
  87    // when this can be turned on or off at runtime, as is the case
  88    // when using the SSE2 registers in DAZ or FTZ mode.
  89    //
  90    static const T m = std::numeric_limits<T>::denorm_min();
  91 #ifdef BOOST_MATH_CHECK_SSE2
  92    return (_mm_getcsr() & (_MM_FLUSH_ZERO_ON | 0x40)) ? tools::min_value<T>() : m;;
  93 #else
  94    return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m;
  95 #endif
  96 }
  97
  98 template <class T>
  99 inline T get_smallest_value(mpl::false_ const&)
 100 {
 101    return tools::min_value<T>();
 102 }
 103
 104 template <class T>
 105 inline T get_smallest_value()
 106 {
 107 #if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310)
 108    return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>());
 109 #else
 110    return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>());
 111 #endif
 112 }
 113
 114 //
 115 // Returns the smallest value that won't generate denorms when
 116 // we calculate the value of the least-significant-bit:
 117 //
 118 template <class T>
 119 T get_min_shift_value();
 120
 121 template <class T>
 122 struct min_shift_initializer
 123 {
 124    struct init
 125    {
 126       init()
 127       {
 128          do_init();
 129       }
 130       static void do_init()
 131       {
 132          get_min_shift_value<T>();
 133       }
 134       void force_instantiate()const{}
 135    };
 136    static const init initializer;
 137    static void force_instantiate()
 138    {
 139       initializer.force_instantiate();
 140    }
 141 };
 142
 143 template <class T>
 144 const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
 145
 146 template <class T>
 147 inline T calc_min_shifted(const mpl::true_&)
 148 {
 149    BOOST_MATH_STD_USING
 150    return ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
 151 }
 152 template <class T>
 153 inline T calc_min_shifted(const mpl::false_&)
 154 {
 155    BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
 156    BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
 157
 158    return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);
 159 }
 160
 161
 162 template <class T>
 163 inline T get_min_shift_value()
 164 {
 165    static const T val = calc_min_shifted<T>(mpl::bool_<!std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());
 166    min_shift_initializer<T>::force_instantiate();
 167
 168    return val;
 169 }
 170
 171 template <class T, class Policy>
 172 T float_next_imp(const T& val, const mpl::true_&, const Policy& pol)
 173 {
 174    BOOST_MATH_STD_USING
 175    int expon;
 176    static const char* function = "float_next<%1%>(%1%)";
 177
 178    int fpclass = (boost::math::fpclassify)(val);
 179
 180    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
 181    {
 182       if(val < 0)
 183          return -tools::max_value<T>();
 184       return policies::raise_domain_error<T>(
 185          function,
 186          "Argument must be finite, but got %1%", val, pol);
 187    }
 188
 189    if(val >= tools::max_value<T>())
 190       return policies::raise_overflow_error<T>(function, 0, pol);
 191
 192    if(val == 0)
 193       return detail::get_smallest_value<T>();
 194
 195    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
 196    {
 197       //
 198       // Special case: if the value of the least significant bit is a denorm, and the result
 199       // would not be a denorm, then shift the input, increment, and shift back.
 200       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
 201       //
 202       return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
 203    }
 204
 205    if(-0.5f == frexp(val, &expon))
 206       --expon; // reduce exponent when val is a power of two, and negative.
 207    T diff = ldexp(T(1), expon - tools::digits<T>());
 208    if(diff == 0)
 209       diff = detail::get_smallest_value<T>();
 210    return val + diff;
 211 } // float_next_imp
 212 //
 213 // Special version for some base other than 2:
 214 //
 215 template <class T, class Policy>
 216 T float_next_imp(const T& val, const mpl::false_&, const Policy& pol)
 217 {
 218    BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
 219    BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
 220
 221    BOOST_MATH_STD_USING
 222    boost::intmax_t expon;
 223    static const char* function = "float_next<%1%>(%1%)";
 224
 225    int fpclass = (boost::math::fpclassify)(val);
 226
 227    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
 228    {
 229       if(val < 0)
 230          return -tools::max_value<T>();
 231       return policies::raise_domain_error<T>(
 232          function,
 233          "Argument must be finite, but got %1%", val, pol);
 234    }
 235
 236    if(val >= tools::max_value<T>())
 237       return policies::raise_overflow_error<T>(function, 0, pol);
 238
 239    if(val == 0)
 240       return detail::get_smallest_value<T>();
 241
 242    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
 243    {
 244       //
 245       // Special case: if the value of the least significant bit is a denorm, and the result
 246       // would not be a denorm, then shift the input, increment, and shift back.
 247       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
 248       //
 249       return scalbn(float_next(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
 250    }
 251
 252    expon = 1 + ilogb(val);
 253    if(-1 == scalbn(val, -expon) * std::numeric_limits<T>::radix)
 254       --expon; // reduce exponent when val is a power of base, and negative.
 255    T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
 256    if(diff == 0)
 257       diff = detail::get_smallest_value<T>();
 258    return val + diff;
 259 } // float_next_imp
 260
 261 } // namespace detail
 262
 263 template <class T, class Policy>
 264 inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
 265 {
 266    typedef typename tools::promote_args<T>::type result_type;
 267    return detail::float_next_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
 268 }
 269
 270 #if 0 //def BOOST_MSVC
 271 //
 272 // We used to use ::_nextafter here, but doing so fails when using
 273 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
 274 // - albeit slower - code instead as at least that gives the correct answer.
 275 //
 276 template <class Policy>
 277 inline double float_next(const double& val, const Policy& pol)
 278 {
 279    static const char* function = "float_next<%1%>(%1%)";
 280
 281    if(!(boost::math::isfinite)(val) && (val > 0))
 282       return policies::raise_domain_error<double>(
 283          function,
 284          "Argument must be finite, but got %1%", val, pol);
 285
 286    if(val >= tools::max_value<double>())
 287       return policies::raise_overflow_error<double>(function, 0, pol);
 288
 289    return ::_nextafter(val, tools::max_value<double>());
 290 }
 291 #endif
 292
 293 template <class T>
 294 inline typename tools::promote_args<T>::type float_next(const T& val)
 295 {
 296    return float_next(val, policies::policy<>());
 297 }
 298
 299 namespace detail{
 300
 301 template <class T, class Policy>
 302 T float_prior_imp(const T& val, const mpl::true_&, const Policy& pol)
 303 {
 304    BOOST_MATH_STD_USING
 305    int expon;
 306    static const char* function = "float_prior<%1%>(%1%)";
 307
 308    int fpclass = (boost::math::fpclassify)(val);
 309
 310    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
 311    {
 312       if(val > 0)
 313          return tools::max_value<T>();
 314       return policies::raise_domain_error<T>(
 315          function,
 316          "Argument must be finite, but got %1%", val, pol);
 317    }
 318
 319    if(val <= -tools::max_value<T>())
 320       return -policies::raise_overflow_error<T>(function, 0, pol);
 321
 322    if(val == 0)
 323       return -detail::get_smallest_value<T>();
 324
 325    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
 326    {
 327       //
 328       // Special case: if the value of the least significant bit is a denorm, and the result
 329       // would not be a denorm, then shift the input, increment, and shift back.
 330       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
 331       //
 332       return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
 333    }
 334
 335    T remain = frexp(val, &expon);
 336    if(remain == 0.5f)
 337       --expon; // when val is a power of two we must reduce the exponent
 338    T diff = ldexp(T(1), expon - tools::digits<T>());
 339    if(diff == 0)
 340       diff = detail::get_smallest_value<T>();
 341    return val - diff;
 342 } // float_prior_imp
 343 //
 344 // Special version for bases other than 2:
 345 //
 346 template <class T, class Policy>
 347 T float_prior_imp(const T& val, const mpl::false_&, const Policy& pol)
 348 {
 349    BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
 350    BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
 351
 352    BOOST_MATH_STD_USING
 353    boost::intmax_t expon;
 354    static const char* function = "float_prior<%1%>(%1%)";
 355
 356    int fpclass = (boost::math::fpclassify)(val);
 357
 358    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
 359    {
 360       if(val > 0)
 361          return tools::max_value<T>();
 362       return policies::raise_domain_error<T>(
 363          function,
 364          "Argument must be finite, but got %1%", val, pol);
 365    }
 366
 367    if(val <= -tools::max_value<T>())
 368       return -policies::raise_overflow_error<T>(function, 0, pol);
 369
 370    if(val == 0)
 371       return -detail::get_smallest_value<T>();
 372
 373    if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
 374    {
 375       //
 376       // Special case: if the value of the least significant bit is a denorm, and the result
 377       // would not be a denorm, then shift the input, increment, and shift back.
 378       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
 379       //
 380       return scalbn(float_prior(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
 381    }
 382
 383    expon = 1 + ilogb(val);
 384    T remain = scalbn(val, -expon);
 385    if(remain * std::numeric_limits<T>::radix == 1)
 386       --expon; // when val is a power of two we must reduce the exponent
 387    T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
 388    if(diff == 0)
 389       diff = detail::get_smallest_value<T>();
 390    return val - diff;
 391 } // float_prior_imp
 392
 393 } // namespace detail
 394
 395 template <class T, class Policy>
 396 inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
 397 {
 398    typedef typename tools::promote_args<T>::type result_type;
 399    return detail::float_prior_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
 400 }
 401
 402 #if 0 //def BOOST_MSVC
 403 //
 404 // We used to use ::_nextafter here, but doing so fails when using
 405 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
 406 // - albeit slower - code instead as at least that gives the correct answer.
 407 //
 408 template <class Policy>
 409 inline double float_prior(const double& val, const Policy& pol)
 410 {
 411    static const char* function = "float_prior<%1%>(%1%)";
 412
 413    if(!(boost::math::isfinite)(val) && (val < 0))
 414       return policies::raise_domain_error<double>(
 415          function,
 416          "Argument must be finite, but got %1%", val, pol);
 417
 418    if(val <= -tools::max_value<double>())
 419       return -policies::raise_overflow_error<double>(function, 0, pol);
 420
 421    return ::_nextafter(val, -tools::max_value<double>());
 422 }
 423 #endif
 424
 425 template <class T>
 426 inline typename tools::promote_args<T>::type float_prior(const T& val)
 427 {
 428    return float_prior(val, policies::policy<>());
 429 }
 430
 431 template <class T, class U, class Policy>
 432 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol)
 433 {
 434    typedef typename tools::promote_args<T, U>::type result_type;
 435    return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol);
 436 }
 437
 438 template <class T, class U>
 439 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
 440 {
 441    return nextafter(val, direction, policies::policy<>());
 442 }
 443
 444 namespace detail{
 445
 446 template <class T, class Policy>
 447 T float_distance_imp(const T& a, const T& b, const mpl::true_&, const Policy& pol)
 448 {
 449    BOOST_MATH_STD_USING
 450    //
 451    // Error handling:
 452    //
 453    static const char* function = "float_distance<%1%>(%1%, %1%)";
 454    if(!(boost::math::isfinite)(a))
 455       return policies::raise_domain_error<T>(
 456          function,
 457          "Argument a must be finite, but got %1%", a, pol);
 458    if(!(boost::math::isfinite)(b))
 459       return policies::raise_domain_error<T>(
 460          function,
 461          "Argument b must be finite, but got %1%", b, pol);
 462    //
 463    // Special cases:
 464    //
 465    if(a > b)
 466       return -float_distance(b, a, pol);
 467    if(a == b)
 468       return T(0);
 469    if(a == 0)
 470       return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
 471    if(b == 0)
 472       return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
 473    if(boost::math::sign(a) != boost::math::sign(b))
 474       return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
 475          + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
 476    //
 477    // By the time we get here, both a and b must have the same sign, we want
 478    // b > a and both postive for the following logic:
 479    //
 480    if(a < 0)
 481       return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
 482
 483    BOOST_ASSERT(a >= 0);
 484    BOOST_ASSERT(b >= a);
 485
 486    int expon;
 487    //
 488    // Note that if a is a denorm then the usual formula fails
 489    // because we actually have fewer than tools::digits<T>()
 490    // significant bits in the representation:
 491    //
 492    (void)frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
 493    T upper = ldexp(T(1), expon);
 494    T result = T(0);
 495    //
 496    // If b is greater than upper, then we *must* split the calculation
 497    // as the size of the ULP changes with each order of magnitude change:
 498    //
 499    if(b > upper)
 500    {
 501       int expon2;
 502       (void)frexp(b, &expon2);
 503       T upper2 = ldexp(T(0.5), expon2);
 504       result = float_distance(upper2, b);
 505       result += (expon2 - expon - 1) * ldexp(T(1), tools::digits<T>() - 1);
 506    }
 507    //
 508    // Use compensated double-double addition to avoid rounding
 509    // errors in the subtraction:
 510    //
 511    expon = tools::digits<T>() - expon;
 512    T mb, x, y, z;
 513    if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
 514    {
 515       //
 516       // Special case - either one end of the range is a denormal, or else the difference is.
 517       // The regular code will fail if we're using the SSE2 registers on Intel and either
 518       // the FTZ or DAZ flags are set.
 519       //
 520       T a2 = ldexp(a, tools::digits<T>());
 521       T b2 = ldexp(b, tools::digits<T>());
 522       mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2);
 523       x = a2 + mb;
 524       z = x - a2;
 525       y = (a2 - (x - z)) + (mb - z);
 526
 527       expon -= tools::digits<T>();
 528    }
 529    else
 530    {
 531       mb = -(std::min)(upper, b);
 532       x = a + mb;
 533       z = x - a;
 534       y = (a - (x - z)) + (mb - z);
 535    }
 536    if(x < 0)
 537    {
 538       x = -x;
 539       y = -y;
 540    }
 541    result += ldexp(x, expon) + ldexp(y, expon);
 542    //
 543    // Result must be an integer:
 544    //
 545    BOOST_ASSERT(result == floor(result));
 546    return result;
 547 } // float_distance_imp
 548 //
 549 // Special versions for bases other than 2:
 550 //
 551 template <class T, class Policy>
 552 T float_distance_imp(const T& a, const T& b, const mpl::false_&, const Policy& pol)
 553 {
 554    BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
 555    BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
 556
 557    BOOST_MATH_STD_USING
 558    //
 559    // Error handling:
 560    //
 561    static const char* function = "float_distance<%1%>(%1%, %1%)";
 562    if(!(boost::math::isfinite)(a))
 563       return policies::raise_domain_error<T>(
 564          function,
 565          "Argument a must be finite, but got %1%", a, pol);
 566    if(!(boost::math::isfinite)(b))
 567       return policies::raise_domain_error<T>(
 568          function,
 569          "Argument b must be finite, but got %1%", b, pol);
 570    //
 571    // Special cases:
 572    //
 573    if(a > b)
 574       return -float_distance(b, a, pol);
 575    if(a == b)
 576       return T(0);
 577    if(a == 0)
 578       return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
 579    if(b == 0)
 580       return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
 581    if(boost::math::sign(a) != boost::math::sign(b))
 582       return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
 583          + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
 584    //
 585    // By the time we get here, both a and b must have the same sign, we want
 586    // b > a and both postive for the following logic:
 587    //
 588    if(a < 0)
 589       return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
 590
 591    BOOST_ASSERT(a >= 0);
 592    BOOST_ASSERT(b >= a);
 593
 594    boost::intmax_t expon;
 595    //
 596    // Note that if a is a denorm then the usual formula fails
 597    // because we actually have fewer than tools::digits<T>()
 598    // significant bits in the representation:
 599    //
 600    expon = 1 + ilogb(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a);
 601    T upper = scalbn(T(1), expon);
 602    T result = T(0);
 603    //
 604    // If b is greater than upper, then we *must* split the calculation
 605    // as the size of the ULP changes with each order of magnitude change:
 606    //
 607    if(b > upper)
 608    {
 609       boost::intmax_t expon2 = 1 + ilogb(b);
 610       T upper2 = scalbn(T(1), expon2 - 1);
 611       result = float_distance(upper2, b);
 612       result += (expon2 - expon - 1) * scalbn(T(1), std::numeric_limits<T>::digits - 1);
 613    }
 614    //
 615    // Use compensated double-double addition to avoid rounding
 616    // errors in the subtraction:
 617    //
 618    expon = std::numeric_limits<T>::digits - expon;
 619    T mb, x, y, z;
 620    if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
 621    {
 622       //
 623       // Special case - either one end of the range is a denormal, or else the difference is.
 624       // The regular code will fail if we're using the SSE2 registers on Intel and either
 625       // the FTZ or DAZ flags are set.
 626       //
 627       T a2 = scalbn(a, std::numeric_limits<T>::digits);
 628       T b2 = scalbn(b, std::numeric_limits<T>::digits);
 629       mb = -(std::min)(T(scalbn(upper, std::numeric_limits<T>::digits)), b2);
 630       x = a2 + mb;
 631       z = x - a2;
 632       y = (a2 - (x - z)) + (mb - z);
 633
 634       expon -= std::numeric_limits<T>::digits;
 635    }
 636    else
 637    {
 638       mb = -(std::min)(upper, b);
 639       x = a + mb;
 640       z = x - a;
 641       y = (a - (x - z)) + (mb - z);
 642    }
 643    if(x < 0)
 644    {
 645       x = -x;
 646       y = -y;
 647    }
 648    result += scalbn(x, expon) + scalbn(y, expon);
 649    //
 650    // Result must be an integer:
 651    //
 652    BOOST_ASSERT(result == floor(result));
 653    return result;
 654 } // float_distance_imp
 655
 656 } // namespace detail
 657
 658 template <class T, class U, class Policy>
 659 inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
 660 {
 661    typedef typename tools::promote_args<T, U>::type result_type;
 662    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()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
 663 }
 664
 665 template <class T, class U>
 666 typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
 667 {
 668    return boost::math::float_distance(a, b, policies::policy<>());
 669 }
 670
 671 namespace detail{
 672
 673 template <class T, class Policy>
 674 T float_advance_imp(T val, int distance, const mpl::true_&, const Policy& pol)
 675 {
 676    BOOST_MATH_STD_USING
 677    //
 678    // Error handling:
 679    //
 680    static const char* function = "float_advance<%1%>(%1%, int)";
 681
 682    int fpclass = (boost::math::fpclassify)(val);
 683
 684    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
 685       return policies::raise_domain_error<T>(
 686          function,
 687          "Argument val must be finite, but got %1%", val, pol);
 688
 689    if(val < 0)
 690       return -float_advance(-val, -distance, pol);
 691    if(distance == 0)
 692       return val;
 693    if(distance == 1)
 694       return float_next(val, pol);
 695    if(distance == -1)
 696       return float_prior(val, pol);
 697
 698    if(fabs(val) < detail::get_min_shift_value<T>())
 699    {
 700       //
 701       // Special case: if the value of the least significant bit is a denorm,
 702       // implement in terms of float_next/float_prior.
 703       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
 704       //
 705       if(distance > 0)
 706       {
 707          do{ val = float_next(val, pol); } while(--distance);
 708       }
 709       else
 710       {
 711          do{ val = float_prior(val, pol); } while(++distance);
 712       }
 713       return val;
 714    }
 715
 716    int expon;
 717    (void)frexp(val, &expon);
 718    T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon);
 719    if(val <= tools::min_value<T>())
 720    {
 721       limit = sign(T(distance)) * tools::min_value<T>();
 722    }
 723    T limit_distance = float_distance(val, limit);
 724    while(fabs(limit_distance) < abs(distance))
 725    {
 726       distance -= itrunc(limit_distance);
 727       val = limit;
 728       if(distance < 0)
 729       {
 730          limit /= 2;
 731          expon--;
 732       }
 733       else
 734       {
 735          limit *= 2;
 736          expon++;
 737       }
 738       limit_distance = float_distance(val, limit);
 739       if(distance && (limit_distance == 0))
 740       {
 741          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);
 742       }
 743    }
 744    if((0.5f == frexp(val, &expon)) && (distance < 0))
 745       --expon;
 746    T diff = 0;
 747    if(val != 0)
 748       diff = distance * ldexp(T(1), expon - tools::digits<T>());
 749    if(diff == 0)
 750       diff = distance * detail::get_smallest_value<T>();
 751    return val += diff;
 752 } // float_advance_imp
 753 //
 754 // Special version for bases other than 2:
 755 //
 756 template <class T, class Policy>
 757 T float_advance_imp(T val, int distance, const mpl::false_&, const Policy& pol)
 758 {
 759    BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
 760    BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
 761
 762    BOOST_MATH_STD_USING
 763    //
 764    // Error handling:
 765    //
 766    static const char* function = "float_advance<%1%>(%1%, int)";
 767
 768    int fpclass = (boost::math::fpclassify)(val);
 769
 770    if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
 771       return policies::raise_domain_error<T>(
 772          function,
 773          "Argument val must be finite, but got %1%", val, pol);
 774
 775    if(val < 0)
 776       return -float_advance(-val, -distance, pol);
 777    if(distance == 0)
 778       return val;
 779    if(distance == 1)
 780       return float_next(val, pol);
 781    if(distance == -1)
 782       return float_prior(val, pol);
 783
 784    if(fabs(val) < detail::get_min_shift_value<T>())
 785    {
 786       //
 787       // Special case: if the value of the least significant bit is a denorm,
 788       // implement in terms of float_next/float_prior.
 789       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
 790       //
 791       if(distance > 0)
 792       {
 793          do{ val = float_next(val, pol); } while(--distance);
 794       }
 795       else
 796       {
 797          do{ val = float_prior(val, pol); } while(++distance);
 798       }
 799       return val;
 800    }
 801
 802    boost::intmax_t expon = 1 + ilogb(val);
 803    T limit = scalbn(T(1), distance < 0 ? expon - 1 : expon);
 804    if(val <= tools::min_value<T>())
 805    {
 806       limit = sign(T(distance)) * tools::min_value<T>();
 807    }
 808    T limit_distance = float_distance(val, limit);
 809    while(fabs(limit_distance) < abs(distance))
 810    {
 811       distance -= itrunc(limit_distance);
 812       val = limit;
 813       if(distance < 0)
 814       {
 815          limit /= std::numeric_limits<T>::radix;
 816          expon--;
 817       }
 818       else
 819       {
 820          limit *= std::numeric_limits<T>::radix;
 821          expon++;
 822       }
 823       limit_distance = float_distance(val, limit);
 824       if(distance && (limit_distance == 0))
 825       {
 826          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);
 827       }
 828    }
 829    /*expon = 1 + ilogb(val);
 830    if((1 == scalbn(val, 1 + expon)) && (distance < 0))
 831       --expon;*/
 832    T diff = 0;
 833    if(val != 0)
 834       diff = distance * scalbn(T(1), expon - std::numeric_limits<T>::digits);
 835    if(diff == 0)
 836       diff = distance * detail::get_smallest_value<T>();
 837    return val += diff;
 838 } // float_advance_imp
 839
 840 } // namespace detail
 841
 842 template <class T, class Policy>
 843 inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
 844 {
 845    typedef typename tools::promote_args<T>::type result_type;
 846    return detail::float_advance_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), distance, mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
 847 }
 848
 849 template <class T>
 850 inline typename tools::promote_args<T>::type float_advance(const T& val, int distance)
 851 {
 852    return boost::math::float_advance(val, distance, policies::policy<>());
 853 }
 854
 855 }} // boost math namespaces
 856
 857 #endif // BOOST_MATH_SPECIAL_NEXT_HPP
 858