diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_log.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_log.hpp index 64cd301e..d1ce8404 100644 --- a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_log.hpp +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_log.hpp @@ -8,28 +8,162 @@ #ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LOG_2024_12_30_HPP #define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LOG_2024_12_30_HPP +#include +#include + #include #include +#if (defined(__GNUC__) && defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128)) +// +// This is the only way we can avoid +// warning: non-standard suffix on floating constant [-Wpedantic] +// when building with -Wall -pedantic. Neither __extension__ +// nor #pragma diagnostic ignored work :( +// +#pragma GCC system_header +#endif + namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { namespace detail { -template -constexpr auto log_impl(T x) -> T +// LCOV_EXCL_START +template +constexpr auto exp_impl(Real x) noexcept -> Real { - // Default to the regular log function. - using std::log; + constexpr int + my_digits_10 + { + ::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits::digits10 + }; - return log(x); + constexpr int + loop_count + { + my_digits_10 < 8 ? 7 + : my_digits_10 < 16 ? 15 + : my_digits_10 < 20 ? 19 + : 33 + }; + + // Scale the argument with a single factor of 2. + // Then square the result upon return. + x /= 2; + + // Perform a simple Taylor series of the exponent function. + Real term { x }; + Real sum { Real { 1 } + term }; + + for (int loop_index { INT8_C(2) }; loop_index < loop_count; ++loop_index) + { + term *= x; + term /= static_cast(loop_index); + sum += term; + } + + // Scale the result. + return sum * sum; } +template +constexpr auto log_impl_pade(Real x) noexcept -> typename std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits::digits < 54), Real>::type +{ + // PadeApproximant[Log[x], {x, 1, {8, 8}}] + // FullSimplify[%] + + const Real + top + { + ((static_cast(-1.0L) + x) * (static_cast(1.0L) + x) * (static_cast(761.0L) + x * (static_cast(28544.0L) + x * (static_cast(209305.0L) + x * (static_cast(423680.0L) + x * (static_cast(209305.0L) + x * (static_cast(28544.0L) + static_cast(761.0L) * x))))))) + }; + + const Real + bot + { + (static_cast(140.0L) * (static_cast(1.0L) + x * (static_cast(64.0L) + x * (static_cast(784.0L) + x * (static_cast(3136.0L) + x * (static_cast(4900.0L) + x * (static_cast(3136.0L) + x * (static_cast(784.0L) + x * (static_cast(64.0L) + x))))))))) + }; + + return top / bot; +} + +template +constexpr auto log_impl_pade(Real x) noexcept -> typename std::enable_if<(!(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits::digits < 54)), Real>::type +{ + // PadeApproximant[Log[x], {x, 1, {16, 16}}] + // FullSimplify[%] + + const Real + top + { + (static_cast(17.0L) * (static_cast(-1.0L) + x) * (static_cast(1.0L) + x) * (static_cast(143327.0L) + x * (static_cast(25160192.0L) + x * (static_cast(1069458527.0L) + x * (static_cast(17931092992.0L) + x * (static_cast(144291009727.0L) + x * (static_cast(613705186816.0L) + x * (static_cast(1446475477311.0L) + x * (static_cast(1923749922816.0L) + x * (static_cast(1446475477311.0L) + x * (static_cast(613705186816.0L) + x * (static_cast(144291009727.0L) + x * (static_cast(17931092992.0L) + x * (static_cast(1069458527.0L) + x * (static_cast(25160192.0L) + static_cast(143327.0L) * x))))))))))))))) + }; + + const Real + bot + { + (static_cast(360360.0L) * (static_cast(1.0L) + x * (static_cast(256.0L) + x * (static_cast(14400.0L) + x * (static_cast(313600.0L) + x * (static_cast(3312400.0L) + x * (static_cast(19079424.0L) + x * (static_cast(64128064.0L) + x * (static_cast(130873600.0L) + x * (static_cast(165636900.0L) + x * (static_cast(130873600.0L) + x * (static_cast(64128064.0L) + x * (static_cast(19079424.0L) + x * (static_cast(3312400.0L) + x * (static_cast(313600.0L) + x * (static_cast(14400.0L) + x * (static_cast(256.0L) + x))))))))))))))))) + }; + + return top / bot; +} + +// N[Log[2], 101] +// 0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754 + +template constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits::digits == 24), FloatingPointType>::type { return static_cast(0.69314718055994530941723212145817656807550013436025525412068L); } +template constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits::digits == 53), FloatingPointType>::type { return static_cast(0.69314718055994530941723212145817656807550013436025525412068L); } +template constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits::digits == 64), FloatingPointType>::type { return static_cast(0.69314718055994530941723212145817656807550013436025525412068L); } +#if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128) +template constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits::digits == 113), FloatingPointType>::type { return static_cast(0.69314718055994530941723212145817656807550013436025525412068Q); } +#else +template constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits::digits == 113), FloatingPointType>::type { return static_cast(0.69314718055994530941723212145817656807550013436025525412068L); } +#endif + +template +constexpr auto log_impl(Real x) noexcept -> Real +{ + int n2 { }; + + Real x2 { ::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::detail::frexp_impl(x, &n2) }; + + if (x2 > static_cast(0.875L)) + { + x2 /= 2; + + ++n2; + } + + Real s { log_impl_pade(x2) }; + + const Real E { exp_impl(s) }; + + // Do one single step of Newton-Raphson iteration. + s = s + ((x2 - E) / E); + + Real xn2 { static_cast(n2) * constant_ln_two() }; + + s += xn2; + + return s; +} +// LCOV_EXCL_STOP + } // namespace detail template constexpr auto log(Real x) -> Real { - return cpp_df_qf_detail::ccmath::detail::log_impl(x); + if (BOOST_MP_IS_CONST_EVALUATED(x)) + { + return detail::log_impl(x); // LCOV_EXCL_LINE + } + else + { + using std::log; + + return log(x); + } } } } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_sqrt.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_sqrt.hpp index 7c9aac5f..3906b989 100644 --- a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_sqrt.hpp +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_sqrt.hpp @@ -52,12 +52,6 @@ constexpr auto sqrt(Real x) noexcept -> Real } } -template -constexpr auto sqrt_constexpr(Real x) -> Real -{ - return detail::sqrt_impl(x); -} - } } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath #endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_SQRT_2023_01_07_HPP diff --git a/include/boost/multiprecision/cpp_double_fp.hpp b/include/boost/multiprecision/cpp_double_fp.hpp index 2c99be49..d58e5203 100644 --- a/include/boost/multiprecision/cpp_double_fp.hpp +++ b/include/boost/multiprecision/cpp_double_fp.hpp @@ -943,26 +943,27 @@ class cpp_double_fp_backend // Use the non-normalized sum of two maximum values, where the lower // value is "shifted" right in the sense of floating-point ldexp. - return cpp_double_fp_backend + return + cpp_double_fp_backend ( arithmetic::two_hilo_sum ( float_type - { - (cpp_df_qf_detail::ccmath::numeric_limits::max)() - * ( - static_cast(1.0F) - - static_cast(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits::epsilon()) - ) - }, + ( + (cpp_df_qf_detail::ccmath::numeric_limits::max)() + * ( + static_cast(1.0F) + - static_cast(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits::epsilon()) + ) + ), float_type - { - cpp_df_qf_detail::ccmath::ldexp - ( - (cpp_df_qf_detail::ccmath::numeric_limits::max)(), - -cpp_df_qf_detail::ccmath::numeric_limits::digits - ) - } + ( + cpp_df_qf_detail::ccmath::ldexp + ( + (cpp_df_qf_detail::ccmath::numeric_limits::max)(), + -cpp_df_qf_detail::ccmath::numeric_limits::digits + ) + ) ) ); } @@ -972,7 +973,8 @@ class cpp_double_fp_backend // Use the non-normalized minimum value, where the lower value // is "shifted" left in the sense of floating-point ldexp. - return cpp_double_fp_backend + return + cpp_double_fp_backend ( float_type ( @@ -987,7 +989,8 @@ class cpp_double_fp_backend static constexpr auto my_value_eps() noexcept -> cpp_double_fp_backend { - return cpp_double_fp_backend + return + cpp_double_fp_backend ( float_type(cpp_df_qf_detail::ccmath::ldexp(float_type { 1 }, int { 3 - my_digits })) ); @@ -1005,12 +1008,40 @@ class cpp_double_fp_backend static constexpr auto my_value_logmax() -> cpp_double_fp_backend { - return float_type(cpp_df_qf_detail::ccmath::log(my_value_max().my_first())); + return + cpp_double_fp_backend + ( + cpp_df_qf_detail::ccmath::log + ( + float_type + ( + (cpp_df_qf_detail::ccmath::numeric_limits::max)() + * ( + static_cast(1.0F) + - static_cast(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits::epsilon()) + ) + ) + ) + ); } static constexpr auto my_value_logmin() -> cpp_double_fp_backend { - return float_type(cpp_df_qf_detail::ccmath::log(my_value_min().my_first())); + return + cpp_double_fp_backend + ( + cpp_df_qf_detail::ccmath::log + ( + float_type + ( + cpp_df_qf_detail::ccmath::ldexp + ( + (cpp_df_qf_detail::ccmath::numeric_limits::min)(), + cpp_df_qf_detail::ccmath::numeric_limits::digits + ) + ) + ) + ); } private: diff --git a/test/test_arithmetic_df.cpp b/test/test_arithmetic_df.cpp index 83bb5348..e21dd540 100644 --- a/test/test_arithmetic_df.cpp +++ b/test/test_arithmetic_df.cpp @@ -30,14 +30,17 @@ auto test_constexpr_ness() noexcept -> void constexpr local_float_type dd_s { sqrt(dd_a) }; constexpr local_float_type sqrt_hundred { sqrt(local_float_type(100)) }; + constexpr local_float_type log_hundred { log(local_float_type(100)) }; static_assert(dd_c > 5, "Error in constexpr multiplication"); static_assert(dd_s > 1, "Error in constexpr square root"); static_assert(sqrt_hundred == 10, "Error in constexpr square root"); + static_assert(log_hundred > local_float_type(4.605L), "Error in constexpr logarithm"); BOOST_CHECK(dd_c > 5); BOOST_CHECK(dd_s > 1); BOOST_CHECK(sqrt_hundred == 10); + BOOST_CHECK(log_hundred > local_float_type(4.605L)); } int main()