Keep the digits of log, log10 and log2 near 1

include/boost/decimal/detail/cmath/impl/log1p_impl.hpp

@@ -9,11 +9,15 @@
 #include <boost/decimal/fwd.hpp>
 #include <boost/decimal/detail/concepts.hpp>
 #include <boost/decimal/detail/cmath/impl/taylor_series_result.hpp>
+#include <boost/decimal/detail/int128.hpp>
+#include <boost/decimal/detail/power_tables.hpp>
+#include <boost/decimal/detail/promotion.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 #include <array>
 #include <cstddef>
 #include <cstdint>
+#include <type_traits>
 #endif
 
 namespace boost {
@@ -200,6 +204,36 @@ constexpr typename log1p_table_imp<b>::d128_fast_coeffs_t log1p_table_imp<b>::d1
 
 using log1p_table = log1p_detail::log1p_table_imp<true>;
 
+// 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 <typename T, typename Array>
+constexpr auto log1p_series_sum(const T z2, const Array& coeffs) noexcept -> T
+{
+    constexpr int digits {std::numeric_limits<T>::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<int>(sig / detail::pow10(static_cast<decltype(sig)>(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<std::size_t>((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>(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>(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>(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>(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>(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>(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 <typename T>
+struct log_scale_t
+{
+    T hi;
+    T lo;
+};
+
+template <typename T, std::enable_if_t<decimal_val_v<T> < 64, bool> = true>
+constexpr auto two_over_ln10() noexcept -> log_scale_t<T>
+{
+    return { T { 8685890, -7 }, T { -3619350, -14 } };
+}
+
+template <typename T, std::enable_if_t<(decimal_val_v<T> >= 64) && (decimal_val_v<T> < 128), bool> = true>
+constexpr auto two_over_ln10() noexcept -> log_scale_t<T>
+{
+    return { T { INT64_C(8685889638065037), -16 }, T { INT64_C(-4469774216216679), -32 } };
+}
+
+template <typename T, std::enable_if_t<decimal_val_v<T> >= 128, bool> = true>
+constexpr auto two_over_ln10() noexcept -> log_scale_t<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 <typename T, std::enable_if_t<decimal_val_v<T> < 64, bool> = true>
+constexpr auto two_over_ln2() noexcept -> log_scale_t<T>
+{
+    return { T { 2885390, -6 }, T { 8177793, -14 } };
+}
+
+template <typename T, std::enable_if_t<(decimal_val_v<T> >= 64) && (decimal_val_v<T> < 128), bool> = true>
+constexpr auto two_over_ln2() noexcept -> log_scale_t<T>
+{
+    return { T { INT64_C(2885390081777927), -15 }, T { INT64_C(-1852801506379962), -31 } };
+}
+
+template <typename T, std::enable_if_t<decimal_val_v<T> >= 128, bool> = true>
+constexpr auto two_over_ln2() noexcept -> log_scale_t<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

include/boost/decimal/detail/cmath/log.hpp

@@ -50,6 +50,11 @@ constexpr auto log_impl(const T x) noexcept
         result = std::numeric_limits<T>::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.

include/boost/decimal/detail/cmath/log10.hpp

@@ -8,6 +8,7 @@
 
 #include <boost/decimal/fwd.hpp> // NOLINT(llvm-include-order)
 #include <boost/decimal/detail/cmath/impl/log_impl.hpp>
+#include <boost/decimal/detail/cmath/log1p.hpp>
 #include <boost/decimal/detail/concepts.hpp>
 #include <boost/decimal/detail/config.hpp>
 #include <boost/decimal/detail/type_traits.hpp>
@@ -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<T>() };
+
+                result = detail::unchecked_fma(parts.wh, k.hi, detail::unchecked_fma(parts.wh, k.lo, parts.s * numbers::log10e_v<T>));
+            }
+            else if (x < one)
             {
                 // Handle reflection.
                 result = -::boost::decimal::log10(one / x);

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 <typename T>
+struct log1p_parts_t
+{
+    T wh;
+    T s;
+};
+
+template <typename T>
+constexpr auto log1p_parts(const T x) noexcept -> log1p_parts_t<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 <typename T>
 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);
         }
     }
 

include/boost/decimal/detail/cmath/log2.hpp

@@ -12,6 +12,7 @@
 #include <boost/decimal/detail/config.hpp>
 #include <boost/decimal/numbers.hpp>
 #include <boost/decimal/detail/cmath/log.hpp>
+#include <boost/decimal/detail/cmath/log1p.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 #include <cmath>
@@ -27,6 +28,17 @@ template <typename T>
 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<T>() };
+
+        return detail::unchecked_fma(parts.wh, k.hi, detail::unchecked_fma(parts.wh, k.lo, parts.s * numbers::log2e_v<T>));
+    }
+
     return ::boost::decimal::log(x) / numbers::ln2_v<T>;
 }
 

test/Jamfile

@@ -140,6 +140,7 @@ run github_issue_1455_downward.cpp
 run github_issue_1459.cpp ;
 run github_issue_1459_toward_zero.cpp
   : : : <pch>off ;
+run github_issue_1467.cpp ;
 
 run link_1.cpp link_2.cpp link_3.cpp ;
 run quick.cpp ;

test/github_issue_1110.cpp

@@ -53,7 +53,7 @@ auto test() -> void
 
         strm << std::setprecision(std::numeric_limits<boost::decimal::decimal128_t>::digits10)<< lgt;
 
-        BOOST_TEST_CSTR_EQ(strm.str().c_str(), "4.4e-33");
+        BOOST_TEST_CSTR_EQ(strm.str().c_str(), "4.342944819032518276511289189166029e-33");
     }
 }
 

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 <boost/decimal.hpp>
+#include <boost/core/lightweight_test.hpp>
+#include <limits>
+
+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 <typename T>
+void check(const T x, const T expected_log, const T expected_log10, const T expected_log2)
+{
+    constexpr T eps {std::numeric_limits<T>::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();
+}