diff --git a/include/boost/decimal/detail/cmath/impl/log1p_impl.hpp b/include/boost/decimal/detail/cmath/impl/log1p_impl.hpp index 161b2295f..72afc03dd 100644 --- a/include/boost/decimal/detail/cmath/impl/log1p_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/log1p_impl.hpp @@ -9,11 +9,15 @@ #include #include #include +#include +#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE #include #include #include +#include #endif namespace boost { @@ -200,6 +204,36 @@ constexpr typename log1p_table_imp::d128_fast_coeffs_t log1p_table_imp::d1 using log1p_table = log1p_detail::log1p_table_imp; +// Horner over the first n terms only. With -log10(z2) >= L, (n + 1) * L >= digits + 2 +// keeps each dropped term less than 0.01 of the last digit of the result. +template +constexpr auto log1p_series_sum(const T z2, const Array& coeffs) noexcept -> T +{ + constexpr int digits {std::numeric_limits::digits10}; + // ceil(1000 * log10(d + 1)) for the leading digit d of z2 + constexpr int log10_next[] {0, 302, 478, 603, 699, 779, 846, 904, 955, 1000}; + + int exp10 {}; + const auto sig {frexp10(z2, &exp10)}; + const int lead {static_cast(sig / detail::pow10(static_cast(digits - 1)))}; + const int bound {1000 * (1 - digits - exp10) - log10_next[lead]}; // 1000 * L + + std::size_t n {coeffs.size()}; + if (bound > 0) + { + const auto need {static_cast((1000 * (digits + 2) + bound - 1) / bound - 1)}; + n = need < 1U ? 1U : (need < n ? need : n); + } + + auto result {coeffs[n - 1U]}; + for (std::size_t i {n - 1U}; i-- > 0U;) + { + result = unchecked_fma(result, z2, coeffs[i]); + } + + return result; +} + // 2*atanh(w) = 2*w + (2/3)*w^3 + (2/5)*w^5 + ... // The first coefficient is exactly 2 for every type, thus the tables leave it out and // the caller adds that term itself. This gives the rest of the series, divided by w^3. @@ -209,37 +243,84 @@ constexpr auto log1p_series_tail(T z2) noexcept; template <> constexpr auto log1p_series_tail(decimal32_t z2) noexcept { - return taylor_series_result(z2, log1p_table::d32_coeffs); + return log1p_series_sum(z2, log1p_table::d32_coeffs); } template <> constexpr auto log1p_series_tail(decimal_fast32_t z2) noexcept { - return taylor_series_result(z2, log1p_table::d32_fast_coeffs); + return log1p_series_sum(z2, log1p_table::d32_fast_coeffs); } template <> constexpr auto log1p_series_tail(decimal64_t z2) noexcept { - return taylor_series_result(z2, log1p_table::d64_coeffs); + return log1p_series_sum(z2, log1p_table::d64_coeffs); } template <> constexpr auto log1p_series_tail(decimal_fast64_t z2) noexcept { - return taylor_series_result(z2, log1p_table::d64_fast_coeffs); + return log1p_series_sum(z2, log1p_table::d64_fast_coeffs); } template <> constexpr auto log1p_series_tail(decimal128_t z2) noexcept { - return taylor_series_result(z2, log1p_table::d128_coeffs); + return log1p_series_sum(z2, log1p_table::d128_coeffs); } template <> constexpr auto log1p_series_tail(decimal_fast128_t z2) noexcept { - return taylor_series_result(z2, log1p_table::d128_fast_coeffs); + return log1p_series_sum(z2, log1p_table::d128_fast_coeffs); +} + +// hi + lo is 2/ln(10) or 2/ln(2) to about two times the digits of the type, with hi +// rounded to the type. log10 and log2 scale the two parts of log1p near 1 with it. +template +struct log_scale_t +{ + T hi; + T lo; +}; + +template < 64, bool> = true> +constexpr auto two_over_ln10() noexcept -> log_scale_t +{ + return { T { 8685890, -7 }, T { -3619350, -14 } }; +} + +template >= 64) && (decimal_val_v < 128), bool> = true> +constexpr auto two_over_ln10() noexcept -> log_scale_t +{ + return { T { INT64_C(8685889638065037), -16 }, T { INT64_C(-4469774216216679), -32 } }; +} + +template >= 128, bool> = true> +constexpr auto two_over_ln10() noexcept -> log_scale_t +{ + return { T { int128::uint128_t { UINT64_C(470863020777972), UINT64_C(4095754970768529350) }, -34 }, + -T { int128::uint128_t { UINT64_C(191964532314735), UINT64_C(3110842849803505067) }, -68 } }; +} + +template < 64, bool> = true> +constexpr auto two_over_ln2() noexcept -> log_scale_t +{ + return { T { 2885390, -6 }, T { 8177793, -14 } }; +} + +template >= 64) && (decimal_val_v < 128), bool> = true> +constexpr auto two_over_ln2() noexcept -> log_scale_t +{ + return { T { INT64_C(2885390081777927), -15 }, T { INT64_C(-1852801506379962), -31 } }; +} + +template >= 128, bool> = true> +constexpr auto two_over_ln2() noexcept -> log_scale_t +{ + return { T { int128::uint128_t { UINT64_C(156417309756587), UINT64_C(14344852839489509192) }, -33 }, + T { int128::uint128_t { UINT64_C(148998268100888), UINT64_C(17076733698656703614) }, -67 } }; } } //namespace detail diff --git a/include/boost/decimal/detail/cmath/log.hpp b/include/boost/decimal/detail/cmath/log.hpp index a8c34928b..88f09286c 100644 --- a/include/boost/decimal/detail/cmath/log.hpp +++ b/include/boost/decimal/detail/cmath/log.hpp @@ -50,6 +50,11 @@ constexpr auto log_impl(const T x) noexcept result = std::numeric_limits::infinity(); } #endif + else if ((x != one) && (x >= T { 5, -1 }) && (x <= T { 15, -1 })) + { + // x - 1 is exact here, and log1p keeps the digits which log10 near 1 cancels. + result = ::boost::decimal::log1p(x - one); + } else if (x < one) { // Handle reflection. diff --git a/include/boost/decimal/detail/cmath/log10.hpp b/include/boost/decimal/detail/cmath/log10.hpp index 468180173..0390b8126 100644 --- a/include/boost/decimal/detail/cmath/log10.hpp +++ b/include/boost/decimal/detail/cmath/log10.hpp @@ -8,6 +8,7 @@ #include // NOLINT(llvm-include-order) #include +#include #include #include #include @@ -76,7 +77,15 @@ constexpr auto log10_impl(const T x) noexcept { constexpr T one { 1 }; - if (x < one) + if ((x >= T { 5, -1 }) && (x <= T { 15, -1 })) + { + // x - 1 is exact here, and log1p keeps the digits which the reduction below cancels. + const auto parts { detail::log1p_parts(x - one) }; + constexpr auto k { detail::two_over_ln10() }; + + result = detail::unchecked_fma(parts.wh, k.hi, detail::unchecked_fma(parts.wh, k.lo, parts.s * numbers::log10e_v)); + } + else if (x < one) { // Handle reflection. result = -::boost::decimal::log10(one / x); diff --git a/include/boost/decimal/detail/cmath/log1p.hpp b/include/boost/decimal/detail/cmath/log1p.hpp index 0ba04be8e..4244d9edc 100644 --- a/include/boost/decimal/detail/cmath/log1p.hpp +++ b/include/boost/decimal/detail/cmath/log1p.hpp @@ -25,6 +25,41 @@ namespace decimal { namespace detail { +// log1p(x) is 2*wh + s for |x| not more than 1/2, before the last rounding. log10 and +// log2 scale the two parts near 1, so that they round only once. +template +struct log1p_parts_t +{ + T wh; + T s; +}; + +template +constexpr auto log1p_parts(const T x) noexcept -> log1p_parts_t +{ + // log1p(x) = 2 * atanh(w), with w = x / (2 + x). For |x| not more than + // 1/2 the value of |w| is not more than 1/3, which the coefficient table + // covers. The first branch takes every larger |x| through log(x + 1). + // Two corrections make the result increase at every argument. The sum + // 2 + x is not exact, and its error makes the quotient fall where x rises, + // thus e holds that error and wh + wl holds the reduction to about two + // times the digits of the type. The first coefficient is exactly 2, thus + // 2*wh stays out of the rounding of the series and only tail rounds. + constexpr T two { 2, 0 }; + + const T d { two + x }; + const T e { x - (d - two) }; // 2 + x == d + e, exact + const T wh { x / d }; + const T r { detail::unchecked_fma(-wh, d, x) - wh * e }; + const T wl { r / d }; + + const T w { wh + wl }; + const T y { w * w }; + const T tail { w * y * detail::log1p_series_tail(y) }; + + return { wh, detail::unchecked_fma(two, wl, tail) }; +} + template constexpr auto log1p_impl(const T x) noexcept BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) @@ -75,27 +110,11 @@ constexpr auto log1p_impl(const T x) noexcept } else { - // log1p(x) = 2 * atanh(w), with w = x / (2 + x). For |x| not more than - // 1/2 the value of |w| is not more than 1/3, which the coefficient table - // covers. The first branch takes every larger |x| through log(x + 1). - // Two corrections make the result increase at every argument. The sum - // 2 + x is not exact, and its error makes the quotient fall where x rises, - // thus e holds that error and wh + wl holds the reduction to about two - // times the digits of the type. The first coefficient is exactly 2, thus - // 2*wh stays out of the rounding of the series and only tail rounds. constexpr T two { 2, 0 }; - const T d { two + x }; - const T e { x - (d - two) }; // 2 + x == d + e, exact - const T wh { x / d }; - const T r { detail::unchecked_fma(-wh, d, x) - wh * e }; - const T wl { r / d }; - - const T w { wh + wl }; - const T y { w * w }; - const T tail { w * y * detail::log1p_series_tail(y) }; + const auto parts { detail::log1p_parts(x) }; - result = detail::unchecked_fma(two, wh, detail::unchecked_fma(two, wl, tail)); + result = detail::unchecked_fma(two, parts.wh, parts.s); } } diff --git a/include/boost/decimal/detail/cmath/log2.hpp b/include/boost/decimal/detail/cmath/log2.hpp index 95fece461..a1b3a8bb5 100644 --- a/include/boost/decimal/detail/cmath/log2.hpp +++ b/include/boost/decimal/detail/cmath/log2.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE #include @@ -27,6 +28,17 @@ template constexpr auto log2_impl(const T x) noexcept BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) { + constexpr T one { 1 }; + + if ((x != one) && (x >= T { 5, -1 }) && (x <= T { 15, -1 })) + { + // As in log10, the two parts of log1p near 1 round only once. + const auto parts { detail::log1p_parts(x - one) }; + constexpr auto k { detail::two_over_ln2() }; + + return detail::unchecked_fma(parts.wh, k.hi, detail::unchecked_fma(parts.wh, k.lo, parts.s * numbers::log2e_v)); + } + return ::boost::decimal::log(x) / numbers::ln2_v; } diff --git a/test/Jamfile b/test/Jamfile index 0c783fbda..9a0f3efa7 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -140,6 +140,7 @@ run github_issue_1455_downward.cpp run github_issue_1459.cpp ; run github_issue_1459_toward_zero.cpp : : : off ; +run github_issue_1467.cpp ; run link_1.cpp link_2.cpp link_3.cpp ; run quick.cpp ; diff --git a/test/github_issue_1110.cpp b/test/github_issue_1110.cpp index f3ee6aaea..ec87818b5 100644 --- a/test/github_issue_1110.cpp +++ b/test/github_issue_1110.cpp @@ -53,7 +53,7 @@ auto test() -> void strm << std::setprecision(std::numeric_limits::digits10)<< lgt; - BOOST_TEST_CSTR_EQ(strm.str().c_str(), "4.4e-33"); + BOOST_TEST_CSTR_EQ(strm.str().c_str(), "4.342944819032518276511289189166029e-33"); } } diff --git a/test/github_issue_1467.cpp b/test/github_issue_1467.cpp new file mode 100644 index 000000000..fb7618324 --- /dev/null +++ b/test/github_issue_1467.cpp @@ -0,0 +1,80 @@ +// Copyright 2026 Shen-Ta Hsieh +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include +#include + +using namespace boost::decimal; +using namespace boost::decimal::literals; + +// The expected values are the MPFR values rounded to the type. Near 1 the old code +// lost most digits, so a relative error of one epsilon still separates the two. +template +void check(const T x, const T expected_log, const T expected_log10, const T expected_log2) +{ + constexpr T eps {std::numeric_limits::epsilon()}; + BOOST_TEST_LE(abs(log(x) - expected_log), eps * abs(expected_log)); + BOOST_TEST_LE(abs(log10(x) - expected_log10), eps * abs(expected_log10)); + BOOST_TEST_LE(abs(log2(x) - expected_log2), eps * abs(expected_log2)); +} + +int main() +{ + check(1.000001_DF, 9.999995e-7_DF, 4.342943e-7_DF, 1.442694e-6_DF); + check(0.9999999_DF, -1.000000e-7_DF, -4.342945e-8_DF, -1.442695e-7_DF); + check(1.000123_DF, 1.229924e-4_DF, 5.341494e-5_DF, 1.774406e-4_DF); + check(0.999377_DF, -6.231941e-4_DF, -2.706498e-4_DF, -8.990791e-4_DF); + check(1.049999_DF, 4.878921e-2_DF, 2.118889e-2_DF, 7.038795e-2_DF); + check(0.9500001_DF, -5.129319e-2_DF, -2.227635e-2_DF, -7.400043e-2_DF); + check(1.499999_DF, 4.054644e-1_DF, 1.760910e-1_DF, 5.849615e-1_DF); + check(0.5000001_DF, -6.931470e-1_DF, -3.010299e-1_DF, -9.999997e-1_DF); + + check(1.000001_DFF, 9.999995e-7_DFF, 4.342943e-7_DFF, 1.442694e-6_DFF); + check(0.9999999_DFF, -1.000000e-7_DFF, -4.342945e-8_DFF, -1.442695e-7_DFF); + check(1.000123_DFF, 1.229924e-4_DFF, 5.341494e-5_DFF, 1.774406e-4_DFF); + check(0.999377_DFF, -6.231941e-4_DFF, -2.706498e-4_DFF, -8.990791e-4_DFF); + check(1.049999_DFF, 4.878921e-2_DFF, 2.118889e-2_DFF, 7.038795e-2_DFF); + check(0.9500001_DFF, -5.129319e-2_DFF, -2.227635e-2_DFF, -7.400043e-2_DFF); + check(1.499999_DFF, 4.054644e-1_DFF, 1.760910e-1_DFF, 5.849615e-1_DFF); + check(0.5000001_DFF, -6.931470e-1_DFF, -3.010299e-1_DFF, -9.999997e-1_DFF); + + check(1.000000000000001_DD, 9.999999999999995e-16_DD, 4.342944819032516e-16_DD, 1.442695040888963e-15_DD); + check(0.9999999999999999_DD, -1.000000000000000e-16_DD, -4.342944819032518e-17_DD, -1.442695040888963e-16_DD); + check(1.000123_DD, 1.229924361202318e-4_DD, 5.341493632285486e-5_DD, 1.774405776575110e-4_DD); + check(0.999377_DD, -6.231941451391355e-4_DD, -2.706497783883408e-4_DD, -8.990791027032677e-4_DD); + check(1.049999_DD, 4.878921178802611e-2_DD, 2.118888545594882e-2_DD, 7.038795389546662e-2_DD); + check(0.9500001_DD, -5.129318912439818e-2_DD, -2.227634899594602e-2_DD, -7.400042958114896e-2_DD); + check(1.499999_DD, 4.054644414412755e-1_DD, 1.760909695259301e-1_DD, 5.849615389241417e-1_DD); + check(0.5000001_DD, -6.931469805599653e-1_DD, -3.010299088050935e-1_DD, -9.999997114610207e-1_DD); + + check(1.000000000000001_DDF, 9.999999999999995e-16_DDF, 4.342944819032516e-16_DDF, 1.442695040888963e-15_DDF); + check(0.9999999999999999_DDF, -1.000000000000000e-16_DDF, -4.342944819032518e-17_DDF, -1.442695040888963e-16_DDF); + check(1.000123_DDF, 1.229924361202318e-4_DDF, 5.341493632285486e-5_DDF, 1.774405776575110e-4_DDF); + check(0.999377_DDF, -6.231941451391355e-4_DDF, -2.706497783883408e-4_DDF, -8.990791027032677e-4_DDF); + check(1.049999_DDF, 4.878921178802611e-2_DDF, 2.118888545594882e-2_DDF, 7.038795389546662e-2_DDF); + check(0.9500001_DDF, -5.129318912439818e-2_DDF, -2.227634899594602e-2_DDF, -7.400042958114896e-2_DDF); + check(1.499999_DDF, 4.054644414412755e-1_DDF, 1.760909695259301e-1_DDF, 5.849615389241417e-1_DDF); + check(0.5000001_DDF, -6.931469805599653e-1_DDF, -3.010299088050935e-1_DDF, -9.999997114610207e-1_DDF); + + check(1.000000000000000000000000000000001_DL, 9.999999999999999999999999999999995e-34_DL, 4.342944819032518276511289189166049e-34_DL, 1.442695040888963407359924681001891e-33_DL); + check(0.9999999999999999999999999999999999_DL, -1.000000000000000000000000000000000e-34_DL, -4.342944819032518276511289189166051e-35_DL, -1.442695040888963407359924681001892e-34_DL); + check(1.000123_DL, 1.229924361202317839697842917749701e-4_DL, 5.341493632285485893172567117842633e-5_DL, 1.774405776575110134764984292573497e-4_DL); + check(0.999377_DL, -6.231941451391354768344471292424107e-4_DL, -2.706497783883407652789107575671064e-4_DL, -8.990791027032676530267330909130347e-4_DL); + check(1.049999_DL, 4.878921178802610708581835504596752e-2_DL, 2.118888545594882490847848964025600e-2_DL, 7.038795389546662003670768661487130e-2_DL); + check(0.9500001_DL, -5.129318912439817885517024157847491e-2_DL, -2.227634899594601834591899834172802e-2_DL, -7.400042958114896379776902251372376e-2_DL); + check(1.499999_DL, 4.054644414412754929903587450939524e-1_DL, 1.760909695259301299856432148730394e-1_DL, 5.849615389241416564376702417853449e-1_DL); + check(0.5000001_DL, -6.931469805599653094145654551915098e-1_DL, -3.010299088050935004518533110905714e-1_DL, -9.999997114610206761042891210845850e-1_DL); + + check(1.000000000000000000000000000000001_DLF, 9.999999999999999999999999999999995e-34_DLF, 4.342944819032518276511289189166049e-34_DLF, 1.442695040888963407359924681001891e-33_DLF); + check(0.9999999999999999999999999999999999_DLF, -1.000000000000000000000000000000000e-34_DLF, -4.342944819032518276511289189166051e-35_DLF, -1.442695040888963407359924681001892e-34_DLF); + check(1.000123_DLF, 1.229924361202317839697842917749701e-4_DLF, 5.341493632285485893172567117842633e-5_DLF, 1.774405776575110134764984292573497e-4_DLF); + check(0.999377_DLF, -6.231941451391354768344471292424107e-4_DLF, -2.706497783883407652789107575671064e-4_DLF, -8.990791027032676530267330909130347e-4_DLF); + check(1.049999_DLF, 4.878921178802610708581835504596752e-2_DLF, 2.118888545594882490847848964025600e-2_DLF, 7.038795389546662003670768661487130e-2_DLF); + check(0.9500001_DLF, -5.129318912439817885517024157847491e-2_DLF, -2.227634899594601834591899834172802e-2_DLF, -7.400042958114896379776902251372376e-2_DLF); + check(1.499999_DLF, 4.054644414412754929903587450939524e-1_DLF, 1.760909695259301299856432148730394e-1_DLF, 5.849615389241416564376702417853449e-1_DLF); + check(0.5000001_DLF, -6.931469805599653094145654551915098e-1_DLF, -3.010299088050935004518533110905714e-1_DLF, -9.999997114610206761042891210845850e-1_DLF); + + return boost::report_errors(); +}