math/ccmath/next.hpp

0001 //  (C) Copyright John Maddock 2008 - 2022.
0002 //  (C) Copyright Matt Borland 2022.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_CCMATH_NEXT_HPP
0008 #define BOOST_MATH_CCMATH_NEXT_HPP
0009
0010 #include <boost/math/ccmath/detail/config.hpp>
0011
0012 #ifdef BOOST_MATH_NO_CCMATH
0013 #error "The header <boost/math/next.hpp> can only be used in C++17 and later."
0014 #endif
0015
0016 #include <stdexcept>
0017 #include <cfloat>
0018 #include <cstdint>
0019 #include <boost/math/policies/policy.hpp>
0020 #include <boost/math/policies/error_handling.hpp>
0021 #include <boost/math/tools/assert.hpp>
0022 #include <boost/math/tools/config.hpp>
0023 #include <boost/math/tools/precision.hpp>
0024 #include <boost/math/tools/traits.hpp>
0025 #include <boost/math/tools/promotion.hpp>
0026 #include <boost/math/ccmath/ilogb.hpp>
0027 #include <boost/math/ccmath/ldexp.hpp>
0028 #include <boost/math/ccmath/scalbln.hpp>
0029 #include <boost/math/ccmath/round.hpp>
0030 #include <boost/math/ccmath/fabs.hpp>
0031 #include <boost/math/ccmath/fpclassify.hpp>
0032 #include <boost/math/ccmath/isfinite.hpp>
0033 #include <boost/math/ccmath/fmod.hpp>
0034
0035 namespace boost::math::ccmath {
0036
0037 namespace detail {
0038
0039 // Forward Declarations
0040 template <typename T, typename result_type = tools::promote_args_t<T>>
0041 constexpr result_type float_prior(const T& val);
0042
0043 template <typename T, typename result_type = tools::promote_args_t<T>>
0044 constexpr result_type float_next(const T& val);
0045
0046 template <typename T>
0047 struct has_hidden_guard_digits;
0048 template <>
0049 struct has_hidden_guard_digits<float> : public std::false_type {};
0050 template <>
0051 struct has_hidden_guard_digits<double> : public std::false_type {};
0052 template <>
0053 struct has_hidden_guard_digits<long double> : public std::false_type {};
0054 #ifdef BOOST_HAS_FLOAT128
0055 template <>
0056 struct has_hidden_guard_digits<__float128> : public std::false_type {};
0057 #endif
0058
0059 template <typename T, bool b>
0060 struct has_hidden_guard_digits_10 : public std::false_type {};
0061 template <typename T>
0062 struct has_hidden_guard_digits_10<T, true> : public std::integral_constant<bool, (std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {};
0063
0064 template <typename T>
0065 struct has_hidden_guard_digits
0066     : public has_hidden_guard_digits_10<T,
0067     std::numeric_limits<T>::is_specialized
0068     && (std::numeric_limits<T>::radix == 10) >
0069 {};
0070
0071 template <typename T>
0072 constexpr T normalize_value(const T& val, const std::false_type&) { return val; }
0073 template <typename T>
0074 constexpr T normalize_value(const T& val, const std::true_type&)
0075 {
0076     static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0077     static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0078
0079     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;
0080     T result = boost::math::ccmath::scalbn(val, shift);
0081     result = boost::math::ccmath::round(result);
0082     return boost::math::ccmath::scalbn(result, -shift);
0083 }
0084
0085 template <typename T>
0086 constexpr T get_smallest_value(const std::true_type&)
0087 {
0088     //
0089     // numeric_limits lies about denorms being present - particularly
0090     // when this can be turned on or off at runtime, as is the case
0091     // when using the SSE2 registers in DAZ or FTZ mode.
0092     //
0093     constexpr T m = std::numeric_limits<T>::denorm_min();
0094     return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m;
0095 }
0096
0097 template <typename T>
0098 constexpr T get_smallest_value(const std::false_type&)
0099 {
0100     return tools::min_value<T>();
0101 }
0102
0103 template <typename T>
0104 constexpr T get_smallest_value()
0105 {
0106     return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized>());
0107 }
0108
0109 template <typename T>
0110 constexpr T calc_min_shifted(const std::true_type&)
0111 {
0112    return boost::math::ccmath::ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
0113 }
0114
0115 template <typename T>
0116 constexpr T calc_min_shifted(const std::false_type&)
0117 {
0118    static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0119    static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0120
0121    return boost::math::ccmath::scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);
0122 }
0123
0124 template <typename T>
0125 constexpr T get_min_shift_value()
0126 {
0127    const T val = calc_min_shifted<T>(std::integral_constant<bool, !std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());
0128    return val;
0129 }
0130
0131 template <typename T, bool b = boost::math::tools::detail::has_backend_type_v<T>>
0132 struct exponent_type
0133 {
0134     using type = int;
0135 };
0136
0137 template <typename T>
0138 struct exponent_type<T, true>
0139 {
0140     using type = typename T::backend_type::exponent_type;
0141 };
0142
0143 template <typename T, bool b = boost::math::tools::detail::has_backend_type_v<T>>
0144 using exponent_type_t = typename exponent_type<T>::type;
0145
0146 template <typename T>
0147 constexpr T float_next_imp(const T& val, const std::true_type&)
0148 {
0149     using exponent_type = exponent_type_t<T>;
0150
0151     exponent_type expon {};
0152
0153     int fpclass = boost::math::ccmath::fpclassify(val);
0154
0155     if (fpclass == FP_NAN)
0156     {
0157         return val;
0158     }
0159     else if (fpclass == FP_INFINITE)
0160     {
0161         return val;
0162     }
0163     else if (val <= -tools::max_value<T>())
0164     {
0165         return val;
0166     }
0167
0168     if (val == 0)
0169     {
0170         return detail::get_smallest_value<T>();
0171     }
0172
0173     if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO)
0174         && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value<T>())
0175         && (val != -tools::min_value<T>()))
0176     {
0177         //
0178         // Special case: if the value of the least significant bit is a denorm, and the result
0179         // would not be a denorm, then shift the input, increment, and shift back.
0180         // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0181         //
0182         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>());
0183     }
0184
0185     if (-0.5f == boost::math::ccmath::frexp(val, &expon))
0186     {
0187         --expon; // reduce exponent when val is a power of two, and negative.
0188     }
0189     T diff = boost::math::ccmath::ldexp(static_cast<T>(1), expon - tools::digits<T>());
0190     if(diff == 0)
0191     {
0192         diff = detail::get_smallest_value<T>();
0193     }
0194     return val + diff;
0195 }
0196
0197 //
0198 // Special version for some base other than 2:
0199 //
0200 template <typename T>
0201 constexpr T float_next_imp(const T& val, const std::false_type&)
0202 {
0203     using exponent_type = exponent_type_t<T>;
0204
0205     static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0206     static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0207
0208     exponent_type expon {};
0209
0210     int fpclass = boost::math::ccmath::fpclassify(val);
0211
0212     if (fpclass == FP_NAN)
0213     {
0214         return val;
0215     }
0216     else if (fpclass == FP_INFINITE)
0217     {
0218         return val;
0219     }
0220     else if (val <= -tools::max_value<T>())
0221     {
0222         return val;
0223     }
0224
0225     if (val == 0)
0226     {
0227         return detail::get_smallest_value<T>();
0228     }
0229
0230     if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO)
0231         && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value<T>())
0232         && (val != -tools::min_value<T>()))
0233     {
0234         //
0235         // Special case: if the value of the least significant bit is a denorm, and the result
0236         // would not be a denorm, then shift the input, increment, and shift back.
0237         // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0238         //
0239         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);
0240     }
0241
0242     expon = 1 + boost::math::ccmath::ilogb(val);
0243     if(-1 == boost::math::ccmath::scalbn(val, -expon) * std::numeric_limits<T>::radix)
0244     {
0245         --expon; // reduce exponent when val is a power of base, and negative.
0246     }
0247
0248     T diff = boost::math::ccmath::scalbn(static_cast<T>(1), expon - std::numeric_limits<T>::digits);
0249     if(diff == 0)
0250     {
0251         diff = detail::get_smallest_value<T>();
0252     }
0253
0254     return val + diff;
0255 }
0256
0257 template <typename T, typename result_type>
0258 constexpr result_type float_next(const T& val)
0259 {
0260     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)>());
0261 }
0262
0263 template <typename T>
0264 constexpr T float_prior_imp(const T& val, const std::true_type&)
0265 {
0266     using exponent_type = exponent_type_t<T>;
0267
0268     exponent_type expon {};
0269
0270     int fpclass = boost::math::ccmath::fpclassify(val);
0271
0272     if (fpclass == FP_NAN)
0273     {
0274         return val;
0275     }
0276     else if (fpclass == FP_INFINITE)
0277     {
0278         return val;
0279     }
0280     else if (val <= -tools::max_value<T>())
0281     {
0282         return val;
0283     }
0284
0285     if (val == 0)
0286     {
0287         return -detail::get_smallest_value<T>();
0288     }
0289
0290     if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO)
0291         && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value<T>())
0292         && (val != tools::min_value<T>()))
0293     {
0294         //
0295         // Special case: if the value of the least significant bit is a denorm, and the result
0296         // would not be a denorm, then shift the input, increment, and shift back.
0297         // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0298         //
0299         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>());
0300     }
0301
0302     if(T remain = boost::math::ccmath::frexp(val, &expon); remain == 0.5f)
0303     {
0304         --expon; // when val is a power of two we must reduce the exponent
0305     }
0306
0307     T diff = boost::math::ccmath::ldexp(static_cast<T>(1), expon - tools::digits<T>());
0308     if(diff == 0)
0309     {
0310         diff = detail::get_smallest_value<T>();
0311     }
0312
0313     return val - diff;
0314 }
0315
0316 //
0317 // Special version for bases other than 2:
0318 //
0319 template <typename T>
0320 constexpr T float_prior_imp(const T& val, const std::false_type&)
0321 {
0322     using exponent_type = exponent_type_t<T>;
0323
0324     static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
0325     static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
0326
0327     exponent_type expon {};
0328
0329     int fpclass = boost::math::ccmath::fpclassify(val);
0330
0331     if (fpclass == FP_NAN)
0332     {
0333         return val;
0334     }
0335     else if (fpclass == FP_INFINITE)
0336     {
0337         return val;
0338     }
0339     else if (val <= -tools::max_value<T>())
0340     {
0341         return val;
0342     }
0343
0344     if (val == 0)
0345     {
0346         return -detail::get_smallest_value<T>();
0347     }
0348
0349     if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO)
0350         && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value<T>())
0351         && (val != tools::min_value<T>()))
0352     {
0353         //
0354         // Special case: if the value of the least significant bit is a denorm, and the result
0355         // would not be a denorm, then shift the input, increment, and shift back.
0356         // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
0357         //
0358         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);
0359     }
0360
0361     expon = 1 + boost::math::ccmath::ilogb(val);
0362
0363     if (T remain = boost::math::ccmath::scalbn(val, -expon); remain * std::numeric_limits<T>::radix == 1)
0364     {
0365         --expon; // when val is a power of two we must reduce the exponent
0366     }
0367
0368     T diff = boost::math::ccmath::scalbn(static_cast<T>(1), expon - std::numeric_limits<T>::digits);
0369     if (diff == 0)
0370     {
0371         diff = detail::get_smallest_value<T>();
0372     }
0373     return val - diff;
0374 } // float_prior_imp
0375
0376 template <typename T, typename result_type>
0377 constexpr result_type float_prior(const T& val)
0378 {
0379     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)>());
0380 }
0381
0382 } // namespace detail
0383
0384 template <typename T, typename U, typename result_type = tools::promote_args_t<T, U>>
0385 constexpr result_type nextafter(const T& val, const U& direction)
0386 {
0387     if (BOOST_MATH_IS_CONSTANT_EVALUATED(val))
0388     {
0389         if (boost::math::ccmath::isnan(val))
0390         {
0391             return val;
0392         }
0393         else if (boost::math::ccmath::isnan(direction))
0394         {
0395             return direction;
0396         }
0397         else if (val < direction)
0398         {
0399             return boost::math::ccmath::detail::float_next(val);
0400         }
0401         else if (val == direction)
0402         {
0403             // IEC 60559 recommends that from is returned whenever from == to. These functions return to instead,
0404             // which makes the behavior around zero consistent: std::nextafter(-0.0, +0.0) returns +0.0 and
0405             // std::nextafter(+0.0, -0.0) returns -0.0.
0406             return direction;
0407         }
0408
0409         return boost::math::ccmath::detail::float_prior(val);
0410     }
0411     else
0412     {
0413         using std::nextafter;
0414         return nextafter(static_cast<result_type>(val), static_cast<result_type>(direction));
0415     }
0416 }
0417
0418 constexpr float nextafterf(float val, float direction)
0419 {
0420     return boost::math::ccmath::nextafter(val, direction);
0421 }
0422
0423 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
0424
0425 constexpr long double nextafterl(long double val, long double direction)
0426 {
0427     return boost::math::ccmath::nextafter(val, direction);
0428 }
0429
0430 template <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>>
0431 constexpr return_type nexttoward(T val, long double direction)
0432 {
0433     if (BOOST_MATH_IS_CONSTANT_EVALUATED(val))
0434     {
0435         return static_cast<return_type>(boost::math::ccmath::nextafter(static_cast<result_type>(val), direction));
0436     }
0437     else
0438     {
0439         using std::nexttoward;
0440         return nexttoward(val, direction);
0441     }
0442 }
0443
0444 constexpr float nexttowardf(float val, long double direction)
0445 {
0446     return boost::math::ccmath::nexttoward(val, direction);
0447 }
0448
0449 constexpr long double nexttowardl(long double val, long double direction)
0450 {
0451     return boost::math::ccmath::nexttoward(val, direction);
0452 }
0453
0454 #endif
0455
0456 } // Namespaces
0457
0458 #endif // BOOST_MATH_SPECIAL_NEXT_HPP