From 54b057bb51f9dee6c780e5ad53deea7ffba84ff2 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sat, 19 Apr 2025 13:35:23 +0100 Subject: [PATCH] Modernise log1p/expm1 to mostly use their std:: equivalents. Add support for types as well. Extend tests for better coverage. --- .../boost/math/special_functions/expm1.hpp | 634 +++++++++--------- .../boost/math/special_functions/log1p.hpp | 115 ++-- include_private/boost/math/tools/test.hpp | 4 +- test/Jamfile.v2 | 2 + test/log1p_expm1_extra_test.cpp | 73 ++ test/log1p_expm1_stdfloat_test.cpp | 87 +++ 6 files changed, 523 insertions(+), 392 deletions(-) create mode 100644 test/log1p_expm1_extra_test.cpp create mode 100644 test/log1p_expm1_stdfloat_test.cpp diff --git a/include/boost/math/special_functions/expm1.hpp b/include/boost/math/special_functions/expm1.hpp index 61809ebcb..651cbb3db 100644 --- a/include/boost/math/special_functions/expm1.hpp +++ b/include/boost/math/special_functions/expm1.hpp @@ -15,6 +15,14 @@ #ifndef BOOST_MATH_HAS_NVRTC +#if defined __has_include +# if __cplusplus > 202002L || _MSVC_LANG > 202002L +# if __has_include () +# include +# endif +# endif +#endif + #include #include #include @@ -37,339 +45,331 @@ #pragma GCC system_header #endif -namespace boost{ namespace math{ +namespace boost { + namespace math { -namespace detail -{ - // Functor expm1_series returns the next term in the Taylor series - // x^k / k! - // each time that operator() is invoked. - // - template - struct expm1_series - { - typedef T result_type; - - BOOST_MATH_GPU_ENABLED expm1_series(T x) - : k(0), m_x(x), m_term(1) {} - - BOOST_MATH_GPU_ENABLED T operator()() - { - ++k; - m_term *= m_x; - m_term /= k; - return m_term; - } - - BOOST_MATH_GPU_ENABLED int count()const - { - return k; - } - - private: - int k; - const T m_x; - T m_term; - expm1_series(const expm1_series&) = delete; - expm1_series& operator=(const expm1_series&) = delete; - }; - -// -// Algorithm expm1 is part of C99, but is not yet provided by many compilers. -// -// This version uses a Taylor series expansion for 0.5 > |x| > epsilon. -// -template -T expm1_imp(T x, const boost::math::integral_constant&, const Policy& pol) -{ - BOOST_MATH_STD_USING - - T a = fabs(x); - if((boost::math::isnan)(a)) - { - return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); - } - if(a > T(0.5f)) - { - if(a >= tools::log_max_value()) + namespace detail { - if(x > 0) - return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, pol); - return -1; - } - return exp(x) - T(1); - } - if(a < tools::epsilon()) - return x; - detail::expm1_series s(x); - boost::math::uintmax_t max_iter = policies::get_max_series_iterations(); + // Functor expm1_series returns the next term in the Taylor series + // x^k / k! + // each time that operator() is invoked. + // + template + struct expm1_series + { + typedef T result_type; - T result = tools::sum_series(s, policies::get_epsilon(), max_iter); + BOOST_MATH_GPU_ENABLED expm1_series(T x) + : k(0), m_x(x), m_term(1) { + } - policies::check_series_iterations("boost::math::expm1<%1%>(%1%)", max_iter, pol); - return result; -} + BOOST_MATH_GPU_ENABLED T operator()() + { + ++k; + m_term *= m_x; + m_term /= k; + return m_term; + } -template -BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant&, const P& pol) -{ - BOOST_MATH_STD_USING + BOOST_MATH_GPU_ENABLED int count()const + { + return k; + } - T a = fabs(x); - if ((boost::math::isnan)(a)) - { - return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); - } - if(a > T(0.5L)) - { - if(a >= tools::log_max_value()) + private: + int k; + const T m_x; + T m_term; + expm1_series(const expm1_series&) = delete; + expm1_series& operator=(const expm1_series&) = delete; + }; + + // + // Algorithm expm1 is part of C99, but is not yet provided by many compilers. + // + // This version uses a Taylor series expansion for 0.5 > |x| > epsilon. + // + template + T expm1_imp(T x, const boost::math::integral_constant&, const Policy& pol) + { + BOOST_MATH_STD_USING + + T a = fabs(x); + if ((boost::math::isnan)(a)) + { + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); + } + if (a > T(0.5f)) + { + if (a >= tools::log_max_value()) + { + if (x > 0) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, pol); + return -1; + } + return exp(x) - T(1); + } + if (a < tools::epsilon()) + return x; + detail::expm1_series s(x); + boost::math::uintmax_t max_iter = policies::get_max_series_iterations(); + + T result = tools::sum_series(s, policies::get_epsilon(), max_iter); + + policies::check_series_iterations("boost::math::expm1<%1%>(%1%)", max_iter, pol); + return result; + } + + template + BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant&, const P& pol) + { + BOOST_MATH_STD_USING + + T a = fabs(x); + if ((boost::math::isnan)(a)) + { + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); + } + if (a > T(0.5L)) + { + if (a >= tools::log_max_value()) + { + if (x > 0) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, pol); + return -1; + } + return exp(x) - T(1); + } + if (a < tools::epsilon()) + return x; + + BOOST_MATH_STATIC const float Y = 0.10281276702880859e1f; + BOOST_MATH_STATIC const T n[] = { static_cast(-0.28127670288085937e-1), static_cast(0.51278186299064534e0), static_cast(-0.6310029069350198e-1), static_cast(0.11638457975729296e-1), static_cast(-0.52143390687521003e-3), static_cast(0.21491399776965688e-4) }; + BOOST_MATH_STATIC const T d[] = { 1, static_cast(-0.45442309511354755e0), static_cast(0.90850389570911714e-1), static_cast(-0.10088963629815502e-1), static_cast(0.63003407478692265e-3), static_cast(-0.17976570003654402e-4) }; + + T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x); + return result; + } + + template + BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant&, const P& pol) + { + BOOST_MATH_STD_USING + + T a = fabs(x); + if ((boost::math::isnan)(a)) + { + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); + } + if (a > T(0.5L)) + { + if (a >= tools::log_max_value()) + { + if (x > 0) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, pol); + return -1; + } + return exp(x) - T(1); + } + if (a < tools::epsilon()) + return x; + + // LCOV_EXCL_START + BOOST_MATH_STATIC const float Y = 0.10281276702880859375e1f; + BOOST_MATH_STATIC const T n[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6) + }; + BOOST_MATH_STATIC const T d[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6) + }; + // LCOV_EXCL_STOP + + T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x); + return result; + } + + template + BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant&, const P& pol) + { + BOOST_MATH_STD_USING + + T a = fabs(x); + if ((boost::math::isnan)(a)) + { + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); + } + if (a > T(0.5L)) + { + if (a >= tools::log_max_value()) + { + if (x > 0) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, pol); + return -1; + } + return exp(x) - T(1); + } + if (a < tools::epsilon()) + return x; + + // LCOV_EXCL_START + static const float Y = 0.10281276702880859375e1f; + static const T n[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10) + }; + static const T d[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12) + }; + // LCOV_EXCL_STOP + + T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x); + return result; + } + + } // namespace detail + + template + BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type expm1(T x, const Policy& /* pol */) { - if(x > 0) - return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, pol); - return -1; + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::precision::type precision_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + typedef boost::math::integral_constant tag_type; + + return policies::checked_narrowing_cast(detail::expm1_imp( + static_cast(x), + tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)"); } - return exp(x) - T(1); - } - if(a < tools::epsilon()) - return x; - BOOST_MATH_STATIC const float Y = 0.10281276702880859e1f; - BOOST_MATH_STATIC const T n[] = { static_cast(-0.28127670288085937e-1), static_cast(0.51278186299064534e0), static_cast(-0.6310029069350198e-1), static_cast(0.11638457975729296e-1), static_cast(-0.52143390687521003e-3), static_cast(0.21491399776965688e-4) }; - BOOST_MATH_STATIC const T d[] = { 1, static_cast(-0.45442309511354755e0), static_cast(0.90850389570911714e-1), static_cast(-0.10088963629815502e-1), static_cast(0.63003407478692265e-3), static_cast(-0.17976570003654402e-4) }; - - T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x); - return result; -} - -template -BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant&, const P& pol) -{ - BOOST_MATH_STD_USING - - T a = fabs(x); - if ((boost::math::isnan)(a)) - { - return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); - } - if(a > T(0.5L)) - { - if(a >= tools::log_max_value()) + // + // Since we now live in a post C++11 world, we can always defer to std::expm1 when appropriate: + // + template + BOOST_MATH_GPU_ENABLED inline float expm1(float x, const Policy&) { - if(x > 0) - return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, pol); - return -1; + BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) + { + if ((boost::math::isnan)(x)) + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); + } + BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) + { + if (x >= tools::log_max_value()) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); + } + return std::expm1(x); } - return exp(x) - T(1); - } - if(a < tools::epsilon()) - return x; - - // LCOV_EXCL_START - BOOST_MATH_STATIC const float Y = 0.10281276702880859375e1f; - BOOST_MATH_STATIC const T n[] = { - BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1), - BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0), - BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1), - BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1), - BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3), - BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4), - BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6) - }; - BOOST_MATH_STATIC const T d[] = { - BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), - BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0), - BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1), - BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1), - BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3), - BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4), - BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6) - }; - // LCOV_EXCL_STOP - - T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x); - return result; -} - -template -BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant&, const P& pol) -{ - BOOST_MATH_STD_USING - - T a = fabs(x); - if ((boost::math::isnan)(a)) - { - return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); - } - if(a > T(0.5L)) - { - if(a >= tools::log_max_value()) +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + template + inline long double expm1(long double x, const Policy&) { - if(x > 0) - return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, pol); - return -1; + BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) + { + if ((boost::math::isnan)(x)) + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); + } + BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) + { + if (x >= tools::log_max_value()) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); + } + return std::expm1(x); } - return exp(x) - T(1); - } - if(a < tools::epsilon()) - return x; - - // LCOV_EXCL_START - static const float Y = 0.10281276702880859375e1f; - static const T n[] = { - BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0), - BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1), - BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4), - BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6), - BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10) - }; - static const T d[] = { - BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), - BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1), - BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2), - BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5), - BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8), - BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10), - BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12) - }; - // LCOV_EXCL_STOP - - T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x); - return result; -} - -} // namespace detail - -template -BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type expm1(T x, const Policy& /* pol */) -{ - typedef typename tools::promote_args::type result_type; - typedef typename policies::evaluation::type value_type; - typedef typename policies::precision::type precision_type; - typedef typename policies::normalise< - Policy, - policies::promote_float, - policies::promote_double, - policies::discrete_quantile<>, - policies::assert_undefined<> >::type forwarding_policy; - - typedef boost::math::integral_constant tag_type; - - return policies::checked_narrowing_cast(detail::expm1_imp( - static_cast(x), - tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)"); -} - -#ifdef expm1 -# ifndef BOOST_HAS_expm1 -# define BOOST_HAS_expm1 -# endif -# undef expm1 #endif + template + BOOST_MATH_GPU_ENABLED inline double expm1(double x, const Policy&) + { + BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) + { + if ((boost::math::isnan)(x)) + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); + } + BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) + { + if (x >= tools::log_max_value()) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); + } + return std::expm1(x); + } -#if defined(BOOST_HAS_EXPM1) && !(defined(__osf__) && defined(__DECCXX_VER)) -# ifdef BOOST_MATH_USE_C99 -template -BOOST_MATH_GPU_ENABLED inline float expm1(float x, const Policy&) -{ - BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) - { - if((boost::math::isnan)(x)) - return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); - } - BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) - { - if (x >= tools::log_max_value()) - return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); - } - return ::expm1f(x); -} -# ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS -template -inline long double expm1(long double x, const Policy&) -{ - BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) - { - if ((boost::math::isnan)(x)) - return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); - } - BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) - { - if (x >= tools::log_max_value()) - return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); - } - return ::expm1l(x); -} -# endif -# else -template -inline float expm1(float x, const Policy&) -{ - BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) - { - if ((boost::math::isnan)(x)) - return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); - } - BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) - { - if (x >= tools::log_max_value()) - return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); - } - return static_cast(::expm1(x)); -} -template -inline typename std::enable_if::type expm1(long double x, const Policy&) -{ - BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) - { - if ((boost::math::isnan)(x)) - return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); - } - BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) - { - if (x >= tools::log_max_value()) - return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); - } - return ::expm1(x); -} -# endif - -template -BOOST_MATH_GPU_ENABLED inline double expm1(double x, const Policy&) -{ - BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) - { - if ((boost::math::isnan)(x)) - return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); - } - BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) - { - if (x >= tools::log_max_value()) - return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); - } - return ::expm1(x); -} + template + BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type expm1(T x) + { + return expm1(x, policies::policy<>()); + } + // + // Specific width floating point types: + // +#ifdef __STDCPP_FLOAT32_T__ + template + BOOST_MATH_GPU_ENABLED inline std::float32_t expm1(std::float32_t x, const Policy& pol) + { + return boost::math::expm1(static_cast(x), pol); + } +#endif +#ifdef __STDCPP_FLOAT64_T__ + template + BOOST_MATH_GPU_ENABLED inline std::float64_t expm1(std::float64_t x, const Policy& pol) + { + return boost::math::expm1(static_cast(x), pol); + } +#endif +#ifdef __STDCPP_FLOAT128_T__ + template + BOOST_MATH_GPU_ENABLED inline std::float128_t expm1(std::float128_t x, const Policy& pol) + { + if constexpr (std::numeric_limits::digits == std::numeric_limits::digits) + { + return boost::math::expm1(static_cast(x), pol); + } + else + { + return boost::math::detail::expm1_imp(x, boost::math::integral_constant(), pol); + } + } #endif - -template -BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type expm1(T x) -{ - return expm1(x, policies::policy<>()); -} - } // namespace math } // namespace boost diff --git a/include/boost/math/special_functions/log1p.hpp b/include/boost/math/special_functions/log1p.hpp index 0366e5632..0fe7ef838 100644 --- a/include/boost/math/special_functions/log1p.hpp +++ b/include/boost/math/special_functions/log1p.hpp @@ -12,6 +12,14 @@ #pragma warning(disable:4702) // Unreachable code (release mode only warning) #endif +#if defined __has_include +# if __cplusplus > 202002L || _MSVC_LANG > 202002L +# if __has_include () +# include +# endif +# endif +#endif + #include #include #include @@ -289,15 +297,6 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type log1p(T x, c detail::log1p_imp(static_cast(x), forwarding_policy(), tag_type()), "boost::math::log1p<%1%>(%1%)"); } -#ifdef log1p -# ifndef BOOST_HAS_LOG1P -# define BOOST_HAS_LOG1P -# endif -# undef log1p -#endif - -#if defined(BOOST_HAS_LOG1P) && !(defined(__osf__) && defined(__DECCXX_VER)) -# ifdef BOOST_MATH_USE_C99 template BOOST_MATH_GPU_ENABLED inline float log1p(float x, const Policy& pol) { @@ -307,7 +306,7 @@ BOOST_MATH_GPU_ENABLED inline float log1p(float x, const Policy& pol) if(x == -1) return -policies::raise_overflow_error( "log1p<%1%>(%1%)", nullptr, pol); - return ::log1pf(x); + return std::log1p(x); } #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS template @@ -319,20 +318,7 @@ BOOST_MATH_GPU_ENABLED inline long double log1p(long double x, const Policy& pol if(x == -1) return -policies::raise_overflow_error( "log1p<%1%>(%1%)", nullptr, pol); - return ::log1pl(x); -} -#endif -#else -template -inline float log1p(float x, const Policy& pol) -{ - if(x < -1) - return policies::raise_domain_error( - "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); - if(x == -1) - return -policies::raise_overflow_error( - "log1p<%1%>(%1%)", nullptr, pol); - return ::log1p(x); + return std::log1p(x); } #endif template @@ -344,56 +330,8 @@ BOOST_MATH_GPU_ENABLED inline double log1p(double x, const Policy& pol) if(x == -1) return -policies::raise_overflow_error( "log1p<%1%>(%1%)", nullptr, pol); - return ::log1p(x); + return std::log1p(x); } -#elif defined(_MSC_VER) && (BOOST_MSVC >= 1400) -// -// You should only enable this branch if you are absolutely sure -// that your compilers optimizer won't mess this code up!! -// Currently tested with VC8 and Intel 9.1. -// -template -inline double log1p(double x, const Policy& pol) -{ - if(x < -1) - return policies::raise_domain_error( - "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); - if(x == -1) - return -policies::raise_overflow_error( - "log1p<%1%>(%1%)", nullptr, pol); - double u = 1+x; - if(u == 1.0) - return x; - else - return ::log(u)*(x/(u-1.0)); -} -template -inline float log1p(float x, const Policy& pol) -{ - return static_cast(boost::math::log1p(static_cast(x), pol)); -} -#ifndef _WIN32_WCE -// -// For some reason this fails to compile under WinCE... -// Needs more investigation. -// -template -inline long double log1p(long double x, const Policy& pol) -{ - if(x < -1) - return policies::raise_domain_error( - "log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol); - if(x == -1) - return -policies::raise_overflow_error( - "log1p<%1%>(%1%)", nullptr, pol); - long double u = 1+x; - if(u == 1.0) - return x; - else - return ::logl(u)*(x/(u-1.0)); -} -#endif -#endif template BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type log1p(T x) @@ -441,6 +379,37 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type log1pmx(T x) return log1pmx(x, policies::policy<>()); } +// +// Specific width floating point types: +// +#ifdef __STDCPP_FLOAT32_T__ +template +BOOST_MATH_GPU_ENABLED inline std::float32_t log1p(std::float32_t x, const Policy& pol) +{ + return boost::math::log1p(static_cast(x), pol); +} +#endif +#ifdef __STDCPP_FLOAT64_T__ +template +BOOST_MATH_GPU_ENABLED inline std::float64_t log1p(std::float64_t x, const Policy& pol) +{ + return boost::math::log1p(static_cast(x), pol); +} +#endif +#ifdef __STDCPP_FLOAT128_T__ +template +BOOST_MATH_GPU_ENABLED inline std::float128_t log1p(std::float128_t x, const Policy& pol) +{ + if constexpr (std::numeric_limits::digits == std::numeric_limits::digits) + { + return boost::math::log1p(static_cast(x), pol); + } + else + { + return boost::math::detail::log1p_imp(x, pol, boost::math::integral_constant()); + } +} +#endif } // namespace math } // namespace boost diff --git a/include_private/boost/math/tools/test.hpp b/include_private/boost/math/tools/test.hpp index 7547ef5be..db7b31616 100644 --- a/include_private/boost/math/tools/test.hpp +++ b/include_private/boost/math/tools/test.hpp @@ -167,7 +167,7 @@ test_result::value_type> test(const A& a, F1 t print_row(row, std::cerr); BOOST_ERROR("Unexpected non-finite result"); } - if(err > 0.5) + if(err > 0.5f) { std::cerr << "CAUTION: Gross error found at entry " << i << ".\n"; std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl; @@ -240,7 +240,7 @@ test_result test_hetero(const A& a, F1 test_func, F2 expect_func) print_row(row, std::cerr); BOOST_ERROR("Unexpected non-finite result"); } - if(err > 0.5) + if(err > 0.5f) { std::cerr << "CAUTION: Gross error found at entry " << i << ".\n"; std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl; diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 47ef67272..2e2a74d6f 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -178,6 +178,8 @@ test-suite special_fun : [ run ccmath_fma_test.cpp /boost/test//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] $(float128_type) ] [ run ccmath_signbit_test.cpp /boost/test//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] $(float128_type) ] [ run log1p_expm1_test.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] + [ run log1p_expm1_extra_test.cpp pch_light /boost/test//boost_unit_test_framework ] + [ run log1p_expm1_stdfloat_test.cpp pch_light /boost/test//boost_unit_test_framework ] [ run powm1_sqrtp1m1_test.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] [ run git_issue_705.cpp /boost/test//boost_unit_test_framework ] [ run git_issue_810.cpp /boost/test//boost_unit_test_framework ] diff --git a/test/log1p_expm1_extra_test.cpp b/test/log1p_expm1_extra_test.cpp new file mode 100644 index 000000000..eb0ff079d --- /dev/null +++ b/test/log1p_expm1_extra_test.cpp @@ -0,0 +1,73 @@ +// Copyright John Maddock 2005. +// Copyright Paul A. Bristow 2010 +// 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 + +#define SC_(x) static_cast(BOOST_STRINGIZE(x)) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#define BOOST_TEST_MAIN +#include +#include +#include +#include +#include "log1p_expm1_test.hpp" +#include + +// +// DESCRIPTION: +// ~~~~~~~~~~~~ +// +// This file tests the functions log1p and expm1. The accuracy tests +// use values generated with NTL::RR at 1000-bit precision +// and our generic versions of these functions. +// +// Note that this tests the generic versions of our functions using +// software emulated types. This ensures these still get tested +// even though we mostly defer to std::expm1 and std::log1p at +// these precisions nowadays +// + +void expected_results() +{ + // + // Define the max and mean errors expected for + // various compilers and platforms. + // + + // + // Catch all cases come last: + // + add_expected_result( + ".*", // compiler + ".*", // stdlib + ".*", // platform + ".*", // test type(s) + ".*", // test data group + ".*", // test function + 4, // Max Peek error + 3); // Max mean error + + // + // Finish off by printing out the compiler/stdlib/platform names, + // we do this to make it easier to mark up expected error rates. + // + std::cout << "Tests run with " << BOOST_COMPILER << ", " + << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl; +} + + +BOOST_AUTO_TEST_CASE( test_main ) +{ + expected_results(); + BOOST_MATH_CONTROL_FP; + test(boost::multiprecision::cpp_bin_float_single(0), "cpp_bin_float_single"); + test(boost::multiprecision::cpp_bin_float_double(0), "cpp_bin_float_double"); + test(boost::multiprecision::cpp_bin_float_double_extended(0), "cpp_bin_float_double_extended"); + test(boost::multiprecision::cpp_bin_float_quad(0), "cpp_bin_float_quad"); +} + diff --git a/test/log1p_expm1_stdfloat_test.cpp b/test/log1p_expm1_stdfloat_test.cpp new file mode 100644 index 000000000..7ad71b4d2 --- /dev/null +++ b/test/log1p_expm1_stdfloat_test.cpp @@ -0,0 +1,87 @@ +// Copyright John Maddock 2005. +// Copyright Paul A. Bristow 2010 +// 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 + +#if defined __has_include +# if __cplusplus > 202002L || _MSVC_LANG > 202002L +# if __has_include () +# include +# endif +# endif +#endif + +#define SC_(x) static_cast(BOOST_MATH_BIG_CONSTANT(T, 128, x)) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#define BOOST_TEST_MAIN +#include +#include +#include +#include +#include +#include "log1p_expm1_test.hpp" + +// +// DESCRIPTION: +// ~~~~~~~~~~~~ +// +// This file tests the functions log1p and expm1. The accuracy tests +// use values generated with NTL::RR at 1000-bit precision +// and our generic versions of these functions. +// +// Note that this tests the generic versions of our functions using +// software emulated types. This ensures these still get tested +// even though we mostly defer to std::expm1 and std::log1p at +// these precisions nowadays +// + +void expected_results() +{ + // + // Define the max and mean errors expected for + // various compilers and platforms. + // + + // + // Catch all cases come last: + // + add_expected_result( + ".*", // compiler + ".*", // stdlib + ".*", // platform + ".*", // test type(s) + ".*", // test data group + ".*", // test function + 4, // Max Peek error + 3); // Max mean error + + // + // Finish off by printing out the compiler/stdlib/platform names, + // we do this to make it easier to mark up expected error rates. + // + std::cout << "Tests run with " << BOOST_COMPILER << ", " + << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl; +} + + +BOOST_AUTO_TEST_CASE( test_main ) +{ + expected_results(); + BOOST_MATH_CONTROL_FP; +#ifdef __STDCPP_FLOAT32_T__ + test(std::float32_t(0), "std::float32_t"); +#endif +#ifdef __STDCPP_FLOAT64_T__ + test(std::float64_t(0), "std::float64_t"); +#endif +#ifdef __STDCPP_FLOAT128_T__ + // As of gcc-13 this does not work as float128_t has no std::math function overloads + //test(std::float128_t(0), "std::float128_t"); +#endif +} +