459 lines · plain
1// (C) Copyright John Maddock 2008 - 2022.2// (C) Copyright Matt Borland 2022.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_CCMATH_NEXT_HPP8#define BOOST_MATH_CCMATH_NEXT_HPP9 10#include <boost/math/ccmath/detail/config.hpp>11 12#ifdef BOOST_MATH_NO_CCMATH13#error "The header <boost/math/next.hpp> can only be used in C++17 and later."14#endif15 16#include <stdexcept>17#include <cfloat>18#include <cstdint>19#include <boost/math/policies/policy.hpp>20#include <boost/math/policies/error_handling.hpp>21#include <boost/math/tools/assert.hpp>22#include <boost/math/tools/config.hpp>23#include <boost/math/tools/precision.hpp>24#include <boost/math/tools/traits.hpp>25#include <boost/math/tools/promotion.hpp>26#include <boost/math/ccmath/ilogb.hpp>27#include <boost/math/ccmath/ldexp.hpp>28#include <boost/math/ccmath/scalbln.hpp>29#include <boost/math/ccmath/round.hpp>30#include <boost/math/ccmath/fabs.hpp>31#include <boost/math/ccmath/fpclassify.hpp>32#include <boost/math/ccmath/isfinite.hpp>33#include <boost/math/ccmath/fmod.hpp>34 35namespace boost::math::ccmath {36 37namespace detail {38 39// Forward Declarations40template <typename T, typename result_type = tools::promote_args_t<T>>41constexpr result_type float_prior(const T& val);42 43template <typename T, typename result_type = tools::promote_args_t<T>>44constexpr result_type float_next(const T& val);45 46template <typename 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#endif58 59template <typename T, bool b>60struct has_hidden_guard_digits_10 : public std::false_type {};61template <typename T>62struct has_hidden_guard_digits_10<T, true> : public std::integral_constant<bool, (std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {};63 64template <typename T>65struct has_hidden_guard_digits 66 : public has_hidden_guard_digits_10<T, 67 std::numeric_limits<T>::is_specialized68 && (std::numeric_limits<T>::radix == 10) >69{};70 71template <typename T>72constexpr T normalize_value(const T& val, const std::false_type&) { return val; }73template <typename T>74constexpr T normalize_value(const T& val, const std::true_type&) 75{76 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");77 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");78 79 std::intmax_t shift = static_cast<std::intmax_t>(std::numeric_limits<T>::digits) - static_cast<std::intmax_t>(boost::math::ccmath::ilogb(val)) - 1;80 T result = boost::math::ccmath::scalbn(val, shift);81 result = boost::math::ccmath::round(result);82 return boost::math::ccmath::scalbn(result, -shift); 83}84 85template <typename T>86constexpr T get_smallest_value(const std::true_type&)87{88 //89 // numeric_limits lies about denorms being present - particularly90 // when this can be turned on or off at runtime, as is the case91 // when using the SSE2 registers in DAZ or FTZ mode.92 //93 constexpr T m = std::numeric_limits<T>::denorm_min();94 return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m;95}96 97template <typename T>98constexpr T get_smallest_value(const std::false_type&)99{100 return tools::min_value<T>();101}102 103template <typename T>104constexpr T get_smallest_value()105{106 return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized>());107}108 109template <typename T>110constexpr T calc_min_shifted(const std::true_type&)111{112 return boost::math::ccmath::ldexp(tools::min_value<T>(), tools::digits<T>() + 1);113}114 115template <typename T>116constexpr T calc_min_shifted(const std::false_type&)117{118 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");119 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");120 121 return boost::math::ccmath::scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);122}123 124template <typename T>125constexpr T get_min_shift_value()126{127 const T val = calc_min_shifted<T>(std::integral_constant<bool, !std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());128 return val;129}130 131template <typename T, bool b = boost::math::tools::detail::has_backend_type_v<T>>132struct exponent_type133{134 using type = int;135};136 137template <typename T>138struct exponent_type<T, true>139{140 using type = typename T::backend_type::exponent_type;141};142 143template <typename T, bool b = boost::math::tools::detail::has_backend_type_v<T>>144using exponent_type_t = typename exponent_type<T>::type;145 146template <typename T>147constexpr T float_next_imp(const T& val, const std::true_type&)148{149 using exponent_type = exponent_type_t<T>;150 151 exponent_type expon {};152 153 int fpclass = boost::math::ccmath::fpclassify(val);154 155 if (fpclass == FP_NAN)156 {157 return val;158 }159 else if (fpclass == FP_INFINITE)160 {161 return val;162 }163 else if (val <= -tools::max_value<T>())164 {165 return val;166 }167 168 if (val == 0)169 {170 return detail::get_smallest_value<T>();171 }172 173 if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) 174 && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value<T>()) 175 && (val != -tools::min_value<T>()))176 {177 //178 // Special case: if the value of the least significant bit is a denorm, and the result179 // would not be a denorm, then shift the input, increment, and shift back.180 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.181 //182 return boost::math::ccmath::ldexp(boost::math::ccmath::detail::float_next(static_cast<T>(boost::math::ccmath::ldexp(val, 2 * tools::digits<T>()))), -2 * tools::digits<T>());183 }184 185 if (-0.5f == boost::math::ccmath::frexp(val, &expon))186 {187 --expon; // reduce exponent when val is a power of two, and negative.188 }189 T diff = boost::math::ccmath::ldexp(static_cast<T>(1), expon - tools::digits<T>());190 if(diff == 0)191 {192 diff = detail::get_smallest_value<T>();193 }194 return val + diff;195}196 197//198// Special version for some base other than 2:199//200template <typename T>201constexpr T float_next_imp(const T& val, const std::false_type&)202{203 using exponent_type = exponent_type_t<T>;204 205 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");206 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");207 208 exponent_type expon {};209 210 int fpclass = boost::math::ccmath::fpclassify(val);211 212 if (fpclass == FP_NAN)213 {214 return val;215 }216 else if (fpclass == FP_INFINITE)217 {218 return val;219 }220 else if (val <= -tools::max_value<T>())221 {222 return val;223 }224 225 if (val == 0)226 {227 return detail::get_smallest_value<T>();228 }229 230 if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) 231 && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value<T>()) 232 && (val != -tools::min_value<T>()))233 {234 //235 // Special case: if the value of the least significant bit is a denorm, and the result236 // would not be a denorm, then shift the input, increment, and shift back.237 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.238 //239 return boost::math::ccmath::scalbn(boost::math::ccmath::detail::float_next(static_cast<T>(boost::math::ccmath::scalbn(val, 2 * std::numeric_limits<T>::digits))), -2 * std::numeric_limits<T>::digits);240 }241 242 expon = 1 + boost::math::ccmath::ilogb(val);243 if(-1 == boost::math::ccmath::scalbn(val, -expon) * std::numeric_limits<T>::radix)244 {245 --expon; // reduce exponent when val is a power of base, and negative.246 }247 248 T diff = boost::math::ccmath::scalbn(static_cast<T>(1), expon - std::numeric_limits<T>::digits);249 if(diff == 0)250 {251 diff = detail::get_smallest_value<T>();252 }253 254 return val + diff;255}256 257template <typename T, typename result_type>258constexpr result_type float_next(const T& val)259{260 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)>());261}262 263template <typename T>264constexpr T float_prior_imp(const T& val, const std::true_type&)265{266 using exponent_type = exponent_type_t<T>;267 268 exponent_type expon {};269 270 int fpclass = boost::math::ccmath::fpclassify(val);271 272 if (fpclass == FP_NAN)273 {274 return val;275 }276 else if (fpclass == FP_INFINITE)277 {278 return val;279 }280 else if (val <= -tools::max_value<T>())281 {282 return val;283 }284 285 if (val == 0)286 {287 return -detail::get_smallest_value<T>();288 }289 290 if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) 291 && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value<T>()) 292 && (val != tools::min_value<T>()))293 {294 //295 // Special case: if the value of the least significant bit is a denorm, and the result296 // would not be a denorm, then shift the input, increment, and shift back.297 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.298 //299 return boost::math::ccmath::ldexp(boost::math::ccmath::detail::float_prior(static_cast<T>(boost::math::ccmath::ldexp(val, 2 * tools::digits<T>()))), -2 * tools::digits<T>());300 }301 302 if(T remain = boost::math::ccmath::frexp(val, &expon); remain == 0.5f)303 {304 --expon; // when val is a power of two we must reduce the exponent305 }306 307 T diff = boost::math::ccmath::ldexp(static_cast<T>(1), expon - tools::digits<T>());308 if(diff == 0)309 {310 diff = detail::get_smallest_value<T>();311 }312 313 return val - diff;314}315 316//317// Special version for bases other than 2:318//319template <typename T>320constexpr T float_prior_imp(const T& val, const std::false_type&)321{322 using exponent_type = exponent_type_t<T>;323 324 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");325 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");326 327 exponent_type expon {};328 329 int fpclass = boost::math::ccmath::fpclassify(val);330 331 if (fpclass == FP_NAN)332 {333 return val;334 }335 else if (fpclass == FP_INFINITE)336 {337 return val;338 }339 else if (val <= -tools::max_value<T>())340 {341 return val;342 }343 344 if (val == 0)345 {346 return -detail::get_smallest_value<T>();347 }348 349 if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) 350 && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value<T>()) 351 && (val != tools::min_value<T>()))352 {353 //354 // Special case: if the value of the least significant bit is a denorm, and the result355 // would not be a denorm, then shift the input, increment, and shift back.356 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.357 //358 return boost::math::ccmath::scalbn(boost::math::ccmath::detail::float_prior(static_cast<T>(boost::math::ccmath::scalbn(val, 2 * std::numeric_limits<T>::digits))), -2 * std::numeric_limits<T>::digits);359 }360 361 expon = 1 + boost::math::ccmath::ilogb(val);362 363 if (T remain = boost::math::ccmath::scalbn(val, -expon); remain * std::numeric_limits<T>::radix == 1)364 {365 --expon; // when val is a power of two we must reduce the exponent366 }367 368 T diff = boost::math::ccmath::scalbn(static_cast<T>(1), expon - std::numeric_limits<T>::digits);369 if (diff == 0)370 {371 diff = detail::get_smallest_value<T>();372 }373 return val - diff;374} // float_prior_imp375 376template <typename T, typename result_type>377constexpr result_type float_prior(const T& val)378{379 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)>());380}381 382} // namespace detail383 384template <typename T, typename U, typename result_type = tools::promote_args_t<T, U>>385constexpr result_type nextafter(const T& val, const U& direction)386{387 if (BOOST_MATH_IS_CONSTANT_EVALUATED(val))388 {389 if (boost::math::ccmath::isnan(val))390 {391 return val;392 }393 else if (boost::math::ccmath::isnan(direction))394 {395 return direction;396 }397 else if (val < direction)398 {399 return boost::math::ccmath::detail::float_next(val);400 }401 else if (val == direction)402 {403 // IEC 60559 recommends that from is returned whenever from == to. These functions return to instead, 404 // which makes the behavior around zero consistent: std::nextafter(-0.0, +0.0) returns +0.0 and 405 // std::nextafter(+0.0, -0.0) returns -0.0.406 return direction;407 }408 409 return boost::math::ccmath::detail::float_prior(val);410 }411 else412 {413 using std::nextafter;414 return nextafter(static_cast<result_type>(val), static_cast<result_type>(direction));415 }416}417 418constexpr float nextafterf(float val, float direction)419{420 return boost::math::ccmath::nextafter(val, direction);421}422 423#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS424 425constexpr long double nextafterl(long double val, long double direction)426{427 return boost::math::ccmath::nextafter(val, direction);428}429 430template <typename T, typename result_type = tools::promote_args_t<T, long double>, typename return_type = std::conditional_t<std::is_integral_v<T>, double, T>>431constexpr return_type nexttoward(T val, long double direction)432{433 if (BOOST_MATH_IS_CONSTANT_EVALUATED(val))434 {435 return static_cast<return_type>(boost::math::ccmath::nextafter(static_cast<result_type>(val), direction));436 }437 else438 {439 using std::nexttoward;440 return nexttoward(val, direction);441 }442}443 444constexpr float nexttowardf(float val, long double direction)445{446 return boost::math::ccmath::nexttoward(val, direction);447}448 449constexpr long double nexttowardl(long double val, long double direction)450{451 return boost::math::ccmath::nexttoward(val, direction);452}453 454#endif455 456} // Namespaces457 458#endif // BOOST_MATH_SPECIAL_NEXT_HPP459