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