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/policies/error_handling.hpp>
  14 #include <boost/math/special_functions/fpclassify.hpp>
  15 #include <boost/math/special_functions/sign.hpp>
  16 #include <boost/math/special_functions/trunc.hpp>
  17
  18 #ifdef BOOST_MSVC
  19 #include <float.h>
  20 #endif
  21
  22 namespace boost{ namespace math{
  23
  24 namespace detail{
  25
  26 template <class T>
  27 inline T get_smallest_value(mpl::true_ const&)
  28 {
  29    //
  30    // numeric_limits lies about denorms being present - particularly
  31    // when this can be turned on or off at runtime, as is the case
  32    // when using the SSE2 registers in DAZ or FTZ mode.
  33    //
  34    static const T m = std::numeric_limits<T>::denorm_min();
  35    return ((tools::min_value<T>() - m) == tools::min_value<T>()) ? tools::min_value<T>() : m;
  36 }
  37
  38 template <class T>
  39 inline T get_smallest_value(mpl::false_ const&)
  40 {
  41    return tools::min_value<T>();
  42 }
  43
  44 template <class T>
  45 inline T get_smallest_value()
  46 {
  47 #if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310)
  48    return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>());
  49 #else
  50    return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>());
  51 #endif
  52 }
  53
  54 //
  55 // Returns the smallest value that won't generate denorms when
  56 // we calculate the value of the least-significant-bit:
  57 //
  58 template <class T>
  59 T get_min_shift_value();
  60
  61 template <class T>
  62 struct min_shift_initializer
  63 {
  64    struct init
  65    {
  66       init()
  67       {
  68          do_init();
  69       }
  70       static void do_init()
  71       {
  72          get_min_shift_value<T>();
  73       }
  74       void force_instantiate()const{}
  75    };
  76    static const init initializer;
  77    static void force_instantiate()
  78    {
  79       initializer.force_instantiate();
  80    }
  81 };
  82
  83 template <class T>
  84 const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
  85
  86
  87 template <class T>
  88 inline T get_min_shift_value()
  89 {
  90    BOOST_MATH_STD_USING
  91    static const T val = ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
  92    min_shift_initializer<T>::force_instantiate();
  93
  94    return val;
  95 }
  96
  97 template <class T, class Policy>
  98 T float_next_imp(const T& val, const Policy& pol)
  99 {
 100    BOOST_MATH_STD_USING
 101    int expon;
 102    static const char* function = "float_next<%1%>(%1%)";
 103
 104    int fpclass = (boost::math::fpclassify)(val);
 105
 106    if((fpclass == FP_NAN) || (fpclass == FP_INFINITE))
 107    {
 108       if(val < 0)
 109          return -tools::max_value<T>();
 110       return policies::raise_domain_error<T>(
 111          function,
 112          "Argument must be finite, but got %1%", val, pol);
 113    }
 114
 115    if(val >= tools::max_value<T>())
 116       return policies::raise_overflow_error<T>(function, 0, pol);
 117
 118    if(val == 0)
 119       return detail::get_smallest_value<T>();
 120
 121    if((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
 122    {
 123       //
 124       // Special case: if the value of the least significant bit is a denorm, and the result
 125       // would not be a denorm, then shift the input, increment, and shift back.
 126       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
 127       //
 128       return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
 129    }
 130
 131    if(-0.5f == frexp(val, &expon))
 132       --expon; // reduce exponent when val is a power of two, and negative.
 133    T diff = ldexp(T(1), expon - tools::digits<T>());
 134    if(diff == 0)
 135       diff = detail::get_smallest_value<T>();
 136    return val + diff;
 137 }
 138
 139 }
 140
 141 template <class T, class Policy>
 142 inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
 143 {
 144    typedef typename tools::promote_args<T>::type result_type;
 145    return detail::float_next_imp(static_cast<result_type>(val), pol);
 146 }
 147
 148 #if 0 //def BOOST_MSVC
 149 //
 150 // We used to use ::_nextafter here, but doing so fails when using
 151 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
 152 // - albeit slower - code instead as at least that gives the correct answer.
 153 //
 154 template <class Policy>
 155 inline double float_next(const double& val, const Policy& pol)
 156 {
 157    static const char* function = "float_next<%1%>(%1%)";
 158
 159    if(!(boost::math::isfinite)(val) && (val > 0))
 160       return policies::raise_domain_error<double>(
 161          function,
 162          "Argument must be finite, but got %1%", val, pol);
 163
 164    if(val >= tools::max_value<double>())
 165       return policies::raise_overflow_error<double>(function, 0, pol);
 166
 167    return ::_nextafter(val, tools::max_value<double>());
 168 }
 169 #endif
 170
 171 template <class T>
 172 inline typename tools::promote_args<T>::type float_next(const T& val)
 173 {
 174    return float_next(val, policies::policy<>());
 175 }
 176
 177 namespace detail{
 178
 179 template <class T, class Policy>
 180 T float_prior_imp(const T& val, const Policy& pol)
 181 {
 182    BOOST_MATH_STD_USING
 183    int expon;
 184    static const char* function = "float_prior<%1%>(%1%)";
 185
 186    int fpclass = (boost::math::fpclassify)(val);
 187
 188    if((fpclass == FP_NAN) || (fpclass == FP_INFINITE))
 189    {
 190       if(val > 0)
 191          return tools::max_value<T>();
 192       return policies::raise_domain_error<T>(
 193          function,
 194          "Argument must be finite, but got %1%", val, pol);
 195    }
 196
 197    if(val <= -tools::max_value<T>())
 198       return -policies::raise_overflow_error<T>(function, 0, pol);
 199
 200    if(val == 0)
 201       return -detail::get_smallest_value<T>();
 202
 203    if((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
 204    {
 205       //
 206       // Special case: if the value of the least significant bit is a denorm, and the result
 207       // would not be a denorm, then shift the input, increment, and shift back.
 208       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
 209       //
 210       return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
 211    }
 212
 213    T remain = frexp(val, &expon);
 214    if(remain == 0.5)
 215       --expon; // when val is a power of two we must reduce the exponent
 216    T diff = ldexp(T(1), expon - tools::digits<T>());
 217    if(diff == 0)
 218       diff = detail::get_smallest_value<T>();
 219    return val - diff;
 220 }
 221
 222 }
 223
 224 template <class T, class Policy>
 225 inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
 226 {
 227    typedef typename tools::promote_args<T>::type result_type;
 228    return detail::float_prior_imp(static_cast<result_type>(val), pol);
 229 }
 230
 231 #if 0 //def BOOST_MSVC
 232 //
 233 // We used to use ::_nextafter here, but doing so fails when using
 234 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
 235 // - albeit slower - code instead as at least that gives the correct answer.
 236 //
 237 template <class Policy>
 238 inline double float_prior(const double& val, const Policy& pol)
 239 {
 240    static const char* function = "float_prior<%1%>(%1%)";
 241
 242    if(!(boost::math::isfinite)(val) && (val < 0))
 243       return policies::raise_domain_error<double>(
 244          function,
 245          "Argument must be finite, but got %1%", val, pol);
 246
 247    if(val <= -tools::max_value<double>())
 248       return -policies::raise_overflow_error<double>(function, 0, pol);
 249
 250    return ::_nextafter(val, -tools::max_value<double>());
 251 }
 252 #endif
 253
 254 template <class T>
 255 inline typename tools::promote_args<T>::type float_prior(const T& val)
 256 {
 257    return float_prior(val, policies::policy<>());
 258 }
 259
 260 template <class T, class U, class Policy>
 261 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol)
 262 {
 263    typedef typename tools::promote_args<T, U>::type result_type;
 264    return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol);
 265 }
 266
 267 template <class T, class U>
 268 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
 269 {
 270    return nextafter(val, direction, policies::policy<>());
 271 }
 272
 273 namespace detail{
 274
 275 template <class T, class Policy>
 276 T float_distance_imp(const T& a, const T& b, const Policy& pol)
 277 {
 278    BOOST_MATH_STD_USING
 279    //
 280    // Error handling:
 281    //
 282    static const char* function = "float_distance<%1%>(%1%, %1%)";
 283    if(!(boost::math::isfinite)(a))
 284       return policies::raise_domain_error<T>(
 285          function,
 286          "Argument a must be finite, but got %1%", a, pol);
 287    if(!(boost::math::isfinite)(b))
 288       return policies::raise_domain_error<T>(
 289          function,
 290          "Argument b must be finite, but got %1%", b, pol);
 291    //
 292    // Special cases:
 293    //
 294    if(a > b)
 295       return -float_distance(b, a, pol);
 296    if(a == b)
 297       return 0;
 298    if(a == 0)
 299       return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
 300    if(b == 0)
 301       return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
 302    if(boost::math::sign(a) != boost::math::sign(b))
 303       return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
 304          + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
 305    //
 306    // By the time we get here, both a and b must have the same sign, we want
 307    // b > a and both postive for the following logic:
 308    //
 309    if(a < 0)
 310       return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
 311
 312    BOOST_ASSERT(a >= 0);
 313    BOOST_ASSERT(b >= a);
 314
 315    int expon;
 316    //
 317    // Note that if a is a denorm then the usual formula fails
 318    // because we actually have fewer than tools::digits<T>()
 319    // significant bits in the representation:
 320    //
 321    frexp(((boost::math::fpclassify)(a) == FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
 322    T upper = ldexp(T(1), expon);
 323    T result = 0;
 324    expon = tools::digits<T>() - expon;
 325    //
 326    // If b is greater than upper, then we *must* split the calculation
 327    // as the size of the ULP changes with each order of magnitude change:
 328    //
 329    if(b > upper)
 330    {
 331       result = float_distance(upper, b);
 332    }
 333    //
 334    // Use compensated double-double addition to avoid rounding
 335    // errors in the subtraction:
 336    //
 337    T mb, x, y, z;
 338    if(((boost::math::fpclassify)(a) == FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
 339    {
 340       //
 341       // Special case - either one end of the range is a denormal, or else the difference is.
 342       // The regular code will fail if we're using the SSE2 registers on Intel and either
 343       // the FTZ or DAZ flags are set.
 344       //
 345       T a2 = ldexp(a, tools::digits<T>());
 346       T b2 = ldexp(b, tools::digits<T>());
 347       mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2);
 348       x = a2 + mb;
 349       z = x - a2;
 350       y = (a2 - (x - z)) + (mb - z);
 351
 352       expon -= tools::digits<T>();
 353    }
 354    else
 355    {
 356       mb = -(std::min)(upper, b);
 357       x = a + mb;
 358       z = x - a;
 359       y = (a - (x - z)) + (mb - z);
 360    }
 361    if(x < 0)
 362    {
 363       x = -x;
 364       y = -y;
 365    }
 366    result += ldexp(x, expon) + ldexp(y, expon);
 367    //
 368    // Result must be an integer:
 369    //
 370    BOOST_ASSERT(result == floor(result));
 371    return result;
 372 }
 373
 374 }
 375
 376 template <class T, class U, class Policy>
 377 inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
 378 {
 379    typedef typename tools::promote_args<T, U>::type result_type;
 380    return detail::float_distance_imp(static_cast<result_type>(a), static_cast<result_type>(b), pol);
 381 }
 382
 383 template <class T, class U>
 384 typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
 385 {
 386    return boost::math::float_distance(a, b, policies::policy<>());
 387 }
 388
 389 namespace detail{
 390
 391 template <class T, class Policy>
 392 T float_advance_imp(T val, int distance, const Policy& pol)
 393 {
 394    BOOST_MATH_STD_USING
 395    //
 396    // Error handling:
 397    //
 398    static const char* function = "float_advance<%1%>(%1%, int)";
 399
 400    int fpclass = (boost::math::fpclassify)(val);
 401
 402    if((fpclass == FP_NAN) || (fpclass == FP_INFINITE))
 403       return policies::raise_domain_error<T>(
 404          function,
 405          "Argument val must be finite, but got %1%", val, pol);
 406
 407    if(val < 0)
 408       return -float_advance(-val, -distance, pol);
 409    if(distance == 0)
 410       return val;
 411    if(distance == 1)
 412       return float_next(val, pol);
 413    if(distance == -1)
 414       return float_prior(val, pol);
 415
 416    if(fabs(val) < detail::get_min_shift_value<T>())
 417    {
 418       //
 419       // Special case: if the value of the least significant bit is a denorm,
 420       // implement in terms of float_next/float_prior.
 421       // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
 422       //
 423       if(distance > 0)
 424       {
 425          do{ val = float_next(val, pol); } while(--distance);
 426       }
 427       else
 428       {
 429          do{ val = float_prior(val, pol); } while(++distance);
 430       }
 431       return val;
 432    }
 433
 434    int expon;
 435    frexp(val, &expon);
 436    T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon);
 437    if(val <= tools::min_value<T>())
 438    {
 439       limit = sign(T(distance)) * tools::min_value<T>();
 440    }
 441    T limit_distance = float_distance(val, limit);
 442    while(fabs(limit_distance) < abs(distance))
 443    {
 444       distance -= itrunc(limit_distance);
 445       val = limit;
 446       if(distance < 0)
 447       {
 448          limit /= 2;
 449          expon--;
 450       }
 451       else
 452       {
 453          limit *= 2;
 454          expon++;
 455       }
 456       limit_distance = float_distance(val, limit);
 457       if(distance && (limit_distance == 0))
 458       {
 459          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);
 460       }
 461    }
 462    if((0.5f == frexp(val, &expon)) && (distance < 0))
 463       --expon;
 464    T diff = 0;
 465    if(val != 0)
 466       diff = distance * ldexp(T(1), expon - tools::digits<T>());
 467    if(diff == 0)
 468       diff = distance * detail::get_smallest_value<T>();
 469    return val += diff;
 470 }
 471
 472 }
 473
 474 template <class T, class Policy>
 475 inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
 476 {
 477    typedef typename tools::promote_args<T>::type result_type;
 478    return detail::float_advance_imp(static_cast<result_type>(val), distance, pol);
 479 }
 480
 481 template <class T>
 482 inline typename tools::promote_args<T>::type float_advance(const T& val, int distance)
 483 {
 484    return boost::math::float_advance(val, distance, policies::policy<>());
 485 }
 486
 487 }} // namespaces
 488
 489 #endif // BOOST_MATH_SPECIAL_NEXT_HPP
 490