From ae0fe3d7514b7a108c8b91b425b22c45336b75ba Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 29 Jun 2022 08:44:54 -0700 Subject: [PATCH] Constexpr next (#789) Implements constexpr: nextafter, nextafterf, nextafterl, nexttoward, nexttowardf, and nexttowardl as described in P0533R9 --- doc/sf/ccmath.qbk | 10 +- include/boost/math/ccmath/next.hpp | 456 ++++++++++++++++++++ include/boost/math/tools/promotion.hpp | 12 + include/boost/math/tools/traits.hpp | 10 + test/Jamfile.v2 | 1 + test/ccmath_next_test.cpp | 62 +++ test/compile_test/ccmath_next_incl_test.cpp | 16 + 7 files changed, 565 insertions(+), 2 deletions(-) create mode 100644 include/boost/math/ccmath/next.hpp create mode 100644 test/ccmath_next_test.cpp create mode 100644 test/compile_test/ccmath_next_incl_test.cpp diff --git a/doc/sf/ccmath.qbk b/doc/sf/ccmath.qbk index 4a2b0b2fbf..6cb6d7df67 100644 --- a/doc/sf/ccmath.qbk +++ b/doc/sf/ccmath.qbk @@ -9,7 +9,7 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) [heading Description] -`Constexpr` implementations of the functionality found in ``. +`Constexpr` implementations of the functionality found in `` and `` [@https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2021/p0533r9.pdf proposed for C++23]. In a `constexpr` context the functions will use an implementation defined in boost. If the context is not `constexpr` the functionality will be directly from the STL implementation of `` used by the compiler. All functions that take an `Integer` type and return a `double` simply cast the `Integer` argument to a `double`. @@ -187,7 +187,13 @@ All of the following functions require C++17 or greater. Requires compiling with fma flag template - inline constepxr Promoted fma(Arithmetic1 x, Arithmetic2 y, Arithmetic3 z) noexcept + inline constexpr Promoted fma(Arithmetic1 x, Arithmetic2 y, Arithmetic3 z) noexcept + + template + constexpr Promoted nextafter(Arithmetic1 from, Arithmetic2 to) + + template + constexpr Promoted nexttoward(T from, long double to) } // Namespaces diff --git a/include/boost/math/ccmath/next.hpp b/include/boost/math/ccmath/next.hpp new file mode 100644 index 0000000000..30045a0a36 --- /dev/null +++ b/include/boost/math/ccmath/next.hpp @@ -0,0 +1,456 @@ +// (C) Copyright John Maddock 2008 - 2022. +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_CCMATH_NEXT_HPP +#define BOOST_MATH_CCMATH_NEXT_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost::math::ccmath { + +namespace detail { + +// Forward Declarations +template > +constexpr result_type float_prior(const T& val); + +template > +constexpr result_type float_next(const T& val); + +template +struct has_hidden_guard_digits; +template <> +struct has_hidden_guard_digits : public std::false_type {}; +template <> +struct has_hidden_guard_digits : public std::false_type {}; +template <> +struct has_hidden_guard_digits : public std::false_type {}; +#ifdef BOOST_HAS_FLOAT128 +template <> +struct has_hidden_guard_digits<__float128> : public std::false_type {}; +#endif + +template +struct has_hidden_guard_digits_10 : public std::false_type {}; +template +struct has_hidden_guard_digits_10 : public std::integral_constant::digits10 != std::numeric_limits::max_digits10)> {}; + +template +struct has_hidden_guard_digits + : public has_hidden_guard_digits_10::is_specialized + && (std::numeric_limits::radix == 10) > +{}; + +template +constexpr T normalize_value(const T& val, const std::false_type&) { return val; } +template +constexpr T normalize_value(const T& val, const std::true_type&) +{ + static_assert(std::numeric_limits::is_specialized, "Type T must be specialized."); + static_assert(std::numeric_limits::radix != 2, "Type T must be specialized."); + + std::intmax_t shift = static_cast(std::numeric_limits::digits) - static_cast(boost::math::ccmath::ilogb(val)) - 1; + T result = boost::math::ccmath::scalbn(val, shift); + result = boost::math::ccmath::round(result); + return boost::math::ccmath::scalbn(result, -shift); +} + +template +constexpr T get_smallest_value(const std::true_type&) +{ + // + // numeric_limits lies about denorms being present - particularly + // when this can be turned on or off at runtime, as is the case + // when using the SSE2 registers in DAZ or FTZ mode. + // + constexpr T m = std::numeric_limits::denorm_min(); + return ((tools::min_value() / 2) == 0) ? tools::min_value() : m; +} + +template +constexpr T get_smallest_value(const std::false_type&) +{ + return tools::min_value(); +} + +template +constexpr T get_smallest_value() +{ + return get_smallest_value(std::integral_constant::is_specialized && (std::numeric_limits::has_denorm == std::denorm_present)>()); +} + +template +constexpr T calc_min_shifted(const std::true_type&) +{ + return boost::math::ccmath::ldexp(tools::min_value(), tools::digits() + 1); +} + +template +constexpr T calc_min_shifted(const std::false_type&) +{ + static_assert(std::numeric_limits::is_specialized, "Type T must be specialized."); + static_assert(std::numeric_limits::radix != 2, "Type T must be specialized."); + + return boost::math::ccmath::scalbn(tools::min_value(), std::numeric_limits::digits + 1); +} + +template +constexpr T get_min_shift_value() +{ + const T val = calc_min_shifted(std::integral_constant::is_specialized || std::numeric_limits::radix == 2>()); + return val; +} + +template > +struct exponent_type +{ + using type = int; +}; + +template +struct exponent_type +{ + using type = typename T::backend_type::exponent_type; +}; + +template > +using exponent_type_t = typename exponent_type::type; + +template +constexpr T float_next_imp(const T& val, const std::true_type&) +{ + using exponent_type = exponent_type_t; + + exponent_type expon {}; + + int fpclass = boost::math::ccmath::fpclassify(val); + + if (fpclass == FP_NAN) + { + return val; + } + else if (fpclass == FP_INFINITE) + { + return val; + } + else if (val <= -tools::max_value()) + { + return val; + } + + if (val == 0) + { + return detail::get_smallest_value(); + } + + if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) + && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value()) + && (val != -tools::min_value())) + { + // + // Special case: if the value of the least significant bit is a denorm, and the result + // would not be a denorm, then shift the input, increment, and shift back. + // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. + // + return boost::math::ccmath::ldexp(boost::math::ccmath::detail::float_next(static_cast(boost::math::ccmath::ldexp(val, 2 * tools::digits()))), -2 * tools::digits()); + } + + if (-0.5f == boost::math::ccmath::frexp(val, &expon)) + { + --expon; // reduce exponent when val is a power of two, and negative. + } + T diff = boost::math::ccmath::ldexp(static_cast(1), expon - tools::digits()); + if(diff == 0) + { + diff = detail::get_smallest_value(); + } + return val + diff; +} + +// +// Special version for some base other than 2: +// +template +constexpr T float_next_imp(const T& val, const std::false_type&) +{ + using exponent_type = exponent_type_t; + + static_assert(std::numeric_limits::is_specialized, "Type T must be specialized."); + static_assert(std::numeric_limits::radix != 2, "Type T must be specialized."); + + exponent_type expon {}; + + int fpclass = boost::math::ccmath::fpclassify(val); + + if (fpclass == FP_NAN) + { + return val; + } + else if (fpclass == FP_INFINITE) + { + return val; + } + else if (val <= -tools::max_value()) + { + return val; + } + + if (val == 0) + { + return detail::get_smallest_value(); + } + + if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) + && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value()) + && (val != -tools::min_value())) + { + // + // Special case: if the value of the least significant bit is a denorm, and the result + // would not be a denorm, then shift the input, increment, and shift back. + // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. + // + return boost::math::ccmath::scalbn(boost::math::ccmath::detail::float_next(static_cast(boost::math::ccmath::scalbn(val, 2 * std::numeric_limits::digits))), -2 * std::numeric_limits::digits); + } + + expon = 1 + boost::math::ccmath::ilogb(val); + if(-1 == boost::math::ccmath::scalbn(val, -expon) * std::numeric_limits::radix) + { + --expon; // reduce exponent when val is a power of base, and negative. + } + + T diff = boost::math::ccmath::scalbn(static_cast(1), expon - std::numeric_limits::digits); + if(diff == 0) + { + diff = detail::get_smallest_value(); + } + + return val + diff; +} + +template +constexpr result_type float_next(const T& val) +{ + return detail::float_next_imp(detail::normalize_value(static_cast(val), typename detail::has_hidden_guard_digits::type()), std::integral_constant::is_specialized || (std::numeric_limits::radix == 2)>()); +} + +template +constexpr T float_prior_imp(const T& val, const std::true_type&) +{ + using exponent_type = exponent_type_t; + + exponent_type expon {}; + + int fpclass = boost::math::ccmath::fpclassify(val); + + if (fpclass == FP_NAN) + { + return val; + } + else if (fpclass == FP_INFINITE) + { + return val; + } + else if (val <= -tools::max_value()) + { + return val; + } + + if (val == 0) + { + return -detail::get_smallest_value(); + } + + if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) + && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value()) + && (val != tools::min_value())) + { + // + // Special case: if the value of the least significant bit is a denorm, and the result + // would not be a denorm, then shift the input, increment, and shift back. + // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. + // + return boost::math::ccmath::ldexp(boost::math::ccmath::detail::float_prior(static_cast(boost::math::ccmath::ldexp(val, 2 * tools::digits()))), -2 * tools::digits()); + } + + if(T remain = boost::math::ccmath::frexp(val, &expon); remain == 0.5f) + { + --expon; // when val is a power of two we must reduce the exponent + } + + T diff = boost::math::ccmath::ldexp(static_cast(1), expon - tools::digits()); + if(diff == 0) + { + diff = detail::get_smallest_value(); + } + + return val - diff; +} + +// +// Special version for bases other than 2: +// +template +constexpr T float_prior_imp(const T& val, const std::false_type&) +{ + using exponent_type = exponent_type_t; + + static_assert(std::numeric_limits::is_specialized, "Type T must be specialized."); + static_assert(std::numeric_limits::radix != 2, "Type T must be specialized."); + + exponent_type expon {}; + + int fpclass = boost::math::ccmath::fpclassify(val); + + if (fpclass == FP_NAN) + { + return val; + } + else if (fpclass == FP_INFINITE) + { + return val; + } + else if (val <= -tools::max_value()) + { + return val; + } + + if (val == 0) + { + return -detail::get_smallest_value(); + } + + if ((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) + && (boost::math::ccmath::fabs(val) < detail::get_min_shift_value()) + && (val != tools::min_value())) + { + // + // Special case: if the value of the least significant bit is a denorm, and the result + // would not be a denorm, then shift the input, increment, and shift back. + // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. + // + return boost::math::ccmath::scalbn(boost::math::ccmath::detail::float_prior(static_cast(boost::math::ccmath::scalbn(val, 2 * std::numeric_limits::digits))), -2 * std::numeric_limits::digits); + } + + expon = 1 + boost::math::ccmath::ilogb(val); + + if (T remain = boost::math::ccmath::scalbn(val, -expon); remain * std::numeric_limits::radix == 1) + { + --expon; // when val is a power of two we must reduce the exponent + } + + T diff = boost::math::ccmath::scalbn(static_cast(1), expon - std::numeric_limits::digits); + if (diff == 0) + { + diff = detail::get_smallest_value(); + } + return val - diff; +} // float_prior_imp + +template +constexpr result_type float_prior(const T& val) +{ + return detail::float_prior_imp(detail::normalize_value(static_cast(val), typename detail::has_hidden_guard_digits::type()), std::integral_constant::is_specialized || (std::numeric_limits::radix == 2)>()); +} + +} // namespace detail + +template > +constexpr result_type nextafter(const T& val, const U& direction) +{ + if (BOOST_MATH_IS_CONSTANT_EVALUATED(val)) + { + if (boost::math::ccmath::isnan(val)) + { + return val; + } + else if (boost::math::ccmath::isnan(direction)) + { + return direction; + } + else if (val < direction) + { + return boost::math::ccmath::detail::float_next(val); + } + else if (val == direction) + { + // IEC 60559 recommends that from is returned whenever from == to. These functions return to instead, + // which makes the behavior around zero consistent: std::nextafter(-0.0, +0.0) returns +0.0 and + // std::nextafter(+0.0, -0.0) returns -0.0. + return direction; + } + + return boost::math::ccmath::detail::float_prior(val); + } + else + { + using std::nextafter; + return nextafter(static_cast(val), static_cast(direction)); + } +} + +constexpr float nextafterf(float val, float direction) +{ + return boost::math::ccmath::nextafter(val, direction); +} + +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + +constexpr long double nextafterl(long double val, long double direction) +{ + return boost::math::ccmath::nextafter(val, direction); +} + +template , typename return_type = std::conditional_t, double, T>> +constexpr return_type nexttoward(T val, long double direction) +{ + if (BOOST_MATH_IS_CONSTANT_EVALUATED(val)) + { + return static_cast(boost::math::ccmath::nextafter(static_cast(val), direction)); + } + else + { + using std::nexttoward; + return nexttoward(val, direction); + } +} + +constexpr float nexttowardf(float val, long double direction) +{ + return boost::math::ccmath::nexttoward(val, direction); +} + +constexpr long double nexttowardl(long double val, long double direction) +{ + return boost::math::ccmath::nexttoward(val, direction); +} + +#endif + +} // Namespaces + +#endif // BOOST_MATH_SPECIAL_NEXT_HPP diff --git a/include/boost/math/tools/promotion.hpp b/include/boost/math/tools/promotion.hpp index 3e51768e6b..72a0755a16 100644 --- a/include/boost/math/tools/promotion.hpp +++ b/include/boost/math/tools/promotion.hpp @@ -63,6 +63,9 @@ namespace boost template <> struct promote_arg { using type = long double; }; template <> struct promote_arg { using type = double; }; + template + using promote_arg_t = typename promote_arg::type; + template struct promote_args_2 { // Promote, if necessary, & pick the wider of the two floating-point types. @@ -108,6 +111,9 @@ namespace boost template <> struct promote_args_2 { using type = long double; }; template <> struct promote_args_2 { using type = long double; }; + template + using promote_args_2_t = typename promote_args_2::type; + template struct promote_args { @@ -135,6 +141,9 @@ namespace boost #endif }; + template + using promote_args_t = typename promote_args::type; + // // This struct is the same as above, but has no static assert on long double usage, // it should be used only on functions that can be implemented for long double @@ -160,6 +169,9 @@ namespace boost >::type; }; + template + using promote_args_permissive_t = typename promote_args_permissive::type; + } // namespace tools } // namespace math } // namespace boost diff --git a/include/boost/math/tools/traits.hpp b/include/boost/math/tools/traits.hpp index cd40f7f1f3..05f6a8632b 100644 --- a/include/boost/math/tools/traits.hpp +++ b/include/boost/math/tools/traits.hpp @@ -53,6 +53,16 @@ BOOST_MATH_HAS_NAMED_TRAIT(has_value_type, value_type) BOOST_MATH_HAS_NAMED_TRAIT(has_policy_type, policy_type) BOOST_MATH_HAS_NAMED_TRAIT(has_backend_type, backend_type) +// C++17-esque helpers +template +constexpr bool has_value_type_v = has_value_type::value; + +template +constexpr bool has_policy_type_v = has_policy_type::value; + +template +constexpr bool has_backend_type_v = has_backend_type::value; + template char cdf(const D& ...); template diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 2708b28096..dba0fcf115 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -159,6 +159,7 @@ test-suite special_fun : [ run ccmath_isless_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] [ run ccmath_islessequal_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] [ run ccmath_isunordered_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] + [ run ccmath_next_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] [ run ccmath_fma_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] [ run log1p_expm1_test.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run powm1_sqrtp1m1_test.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] diff --git a/test/ccmath_next_test.cpp b/test/ccmath_next_test.cpp new file mode 100644 index 0000000000..881155f27e --- /dev/null +++ b/test/ccmath_next_test.cpp @@ -0,0 +1,62 @@ +// (C) Copyright John Maddock 2008 - 2022. +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include +#include +#include "math_unit_test.hpp" + +#if !defined(BOOST_MATH_NO_CONSTEXPR_DETECTION) && !defined(BOOST_MATH_USING_BUILTIN_CONSTANT_P) +template +void test_next() +{ + // NaN handling + static_assert(boost::math::ccmath::isnan(boost::math::ccmath::nextafter(std::numeric_limits::quiet_NaN(), T(0)))); + static_assert(boost::math::ccmath::isnan(boost::math::ccmath::nextafter(T(0), std::numeric_limits::quiet_NaN()))); + + // Handling of 0 + static_assert(boost::math::ccmath::nextafter(T(-0.0), T(0.0)) == T(0.0)); + static_assert(boost::math::ccmath::nextafter(T(0.0), T(-0.0)) == T(-0.0)); + + // val = 1 + constexpr T test_1 = boost::math::ccmath::nextafter(T(1), T(1.5)); + static_assert(test_1 < 1 + 2*std::numeric_limits::epsilon()); + static_assert(test_1 > 1 - 2*std::numeric_limits::epsilon()); + + constexpr T test_1_toward = boost::math::ccmath::nexttoward(T(1), T(1.5)); + + // For T is long double nextafter is the same as nexttoward + // For T is not long double the answer will be either greater or equal when from > to depending on loss of precision + static_assert(test_1 >= test_1_toward); + + // Compare to existing implementation + // test_1 has already passed through static_asserts so we know it was calculated at compile time + // rather than farming out to std at run time. + const T existing_test_1 = boost::math::nextafter(T(1), T(1.5)); + CHECK_EQUAL(test_1, existing_test_1); +} + +int main(void) +{ + test_next(); + test_next(); + + #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_next(); + #endif + + return boost::math::test::report_errors(); +} +#else +int main(void) +{ + return 0; +} +#endif diff --git a/test/compile_test/ccmath_next_incl_test.cpp b/test/compile_test/ccmath_next_incl_test.cpp new file mode 100644 index 0000000000..d9ee72b227 --- /dev/null +++ b/test/compile_test/ccmath_next_incl_test.cpp @@ -0,0 +1,16 @@ +// (C) Copyright Matt Borland 2022. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include "test_compile_result.hpp" + +void compile_and_link_test() +{ + check_result(boost::math::ccmath::nextafter(1.0F, 1.05F)); + check_result(boost::math::ccmath::nextafter(1.0, 1.0)); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + check_result(boost::math::ccmath::nexttoward(1.0L, 1.0L)); +#endif +}