diff --git a/libc/src/math/generic/log1p.cpp b/libc/src/math/generic/log1p.cpp index 43eb8a924aef4..b9c58b843a240 100644 --- a/libc/src/math/generic/log1p.cpp +++ b/libc/src/math/generic/log1p.cpp @@ -822,8 +822,8 @@ constexpr Float128 BIG_COEFFS[4]{ {Sign::NEG, -128, 0x80000000'00000000'00000000'00000000_u128}, }; -LIBC_INLINE double log1p_accurate(int e_x, int index, - fputil::DoubleDouble m_x) { +[[maybe_unused]] LIBC_INLINE double log1p_accurate(int e_x, int index, + fputil::DoubleDouble m_x) { Float128 e_x_f128(static_cast(e_x)); Float128 sum = fputil::quick_mul(LOG_2, e_x_f128); sum = fputil::quick_add(sum, LOG_R1[index]); @@ -882,7 +882,6 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) { constexpr int EXP_BIAS = FPBits_t::EXP_BIAS; constexpr int FRACTION_LEN = FPBits_t::FRACTION_LEN; - constexpr uint64_t FRACTION_MASK = FPBits_t::FRACTION_MASK; FPBits_t xbits(x); uint64_t x_u = xbits.uintval(); @@ -954,12 +953,12 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) { // |x_dd.lo| < ulp(x_dd.hi) FPBits_t xhi_bits(x_dd.hi); + uint64_t xhi_frac = xhi_bits.get_mantissa(); x_u = xhi_bits.uintval(); // Range reduction: // Find k such that |x_hi - k * 2^-7| <= 2^-8. - int idx = - static_cast(((x_u & FRACTION_MASK) + (1ULL << (FRACTION_LEN - 8))) >> - (FRACTION_LEN - 7)); + int idx = static_cast((xhi_frac + (1ULL << (FRACTION_LEN - 8))) >> + (FRACTION_LEN - 7)); int x_e = xhi_bits.get_exponent() + (idx >> 7); double e_x = static_cast(x_e); @@ -974,17 +973,21 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) { constexpr double ERR_HI[2] = {0x1.0p-85, 0.0}; double err_hi = ERR_HI[hi == 0.0]; - // Scaling factior = 2^(-xh_bits.get_exponent()) - uint64_t s_u = (static_cast(EXP_BIAS) << (FRACTION_LEN + 1)) - - (x_u & FPBits_t::EXP_MASK); - // When the exponent of x is 2^1023, its inverse, 2^(-1023), is subnormal. - const double EXPONENT_CORRECTION[2] = {0.0, 0x1.0p-1023}; - double scaling = FPBits_t(s_u).get_val() + EXPONENT_CORRECTION[s_u == 0]; + // Scale x_dd by 2^(-xh_bits.get_exponent()). + int64_t s_u = static_cast(x_u & FPBits_t::EXP_MASK) - + (static_cast(EXP_BIAS) << FRACTION_LEN); // Normalize arguments: // 1 <= m_dd.hi < 2 // |m_dd.lo| < 2^-52. // This is exact. - fputil::DoubleDouble m_dd{scaling * x_dd.lo, scaling * x_dd.hi}; + uint64_t m_hi = FPBits_t::one().uintval() | xhi_frac; + + uint64_t m_lo = + FPBits_t(x_dd.lo).abs().get_val() > x_dd.hi * 0x1.0p-127 + ? static_cast(cpp::bit_cast(x_dd.lo) - s_u) + : 0; + + fputil::DoubleDouble m_dd{FPBits_t(m_lo).get_val(), FPBits_t(m_hi).get_val()}; // Perform range reduction: // r * m - 1 = r * (m_dd.hi + m_dd.lo) - 1 diff --git a/libc/test/src/math/smoke/log1p_test.cpp b/libc/test/src/math/smoke/log1p_test.cpp index eba65f56df739..b98c0f26a8bca 100644 --- a/libc/test/src/math/smoke/log1p_test.cpp +++ b/libc/test/src/math/smoke/log1p_test.cpp @@ -13,8 +13,6 @@ #include "test/UnitTest/FPMatcher.h" #include "test/UnitTest/Test.h" -#include - using LlvmLibcLog1pTest = LIBC_NAMESPACE::testing::FPTest; TEST_F(LlvmLibcLog1pTest, SpecialNumbers) { @@ -26,6 +24,9 @@ TEST_F(LlvmLibcLog1pTest, SpecialNumbers) { EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::log1p(-0.0)); EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, LIBC_NAMESPACE::log1p(-1.0), FE_DIVBYZERO); + + EXPECT_FP_EQ(0x1.62c829bf8fd9dp9, + LIBC_NAMESPACE::log1p(0x1.9b536cac3a09dp1023)); } #ifdef LIBC_TEST_FTZ_DAZ @@ -36,18 +37,24 @@ TEST_F(LlvmLibcLog1pTest, FTZMode) { ModifyMXCSR mxcsr(FTZ); EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::log1p(min_denormal)); + EXPECT_FP_EQ(0x1.62c829bf8fd9dp9, + LIBC_NAMESPACE::log1p(0x1.9b536cac3a09dp1023)); } TEST_F(LlvmLibcLog1pTest, DAZMode) { ModifyMXCSR mxcsr(DAZ); EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::log1p(min_denormal)); + EXPECT_FP_EQ(0x1.62c829bf8fd9dp9, + LIBC_NAMESPACE::log1p(0x1.9b536cac3a09dp1023)); } TEST_F(LlvmLibcLog1pTest, FTZDAZMode) { ModifyMXCSR mxcsr(FTZ | DAZ); EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::log1p(min_denormal)); + EXPECT_FP_EQ(0x1.62c829bf8fd9dp9, + LIBC_NAMESPACE::log1p(0x1.9b536cac3a09dp1023)); } #endif