[libc][math] Fix log1p SEGV with large inputs when FTZ/DAZ flags are set.

[lntue]

No description provided.

[llvmbot]

@llvm/pr-subscribers-libc

Author: None (lntue)

Changes

---

Full diff: https://github.com/llvm/llvm-project/pull/115541.diff

2 Files Affected:

• (modified) libc/src/math/generic/log1p.cpp (+7-7)
• (modified) libc/test/src/math/smoke/log1p_test.cpp (+5-2)

diff --git a/libc/src/math/generic/log1p.cpp b/libc/src/math/generic/log1p.cpp
index 43eb8a924aef47..ec1f2bc00fd81f 100644
--- a/libc/src/math/generic/log1p.cpp
+++ b/libc/src/math/generic/log1p.cpp
@@ -822,7 +822,7 @@ constexpr Float128 BIG_COEFFS[4]{
     {Sign::NEG, -128, 0x80000000'00000000'00000000'00000000_u128},
 };
 
-LIBC_INLINE double log1p_accurate(int e_x, int index,
+[[maybe_unused]] LIBC_INLINE double log1p_accurate(int e_x, int index,
                                   fputil::DoubleDouble m_x) {
   Float128 e_x_f128(static_cast<float>(e_x));
   Float128 sum = fputil::quick_mul(LOG_2, e_x_f128);
@@ -975,16 +975,16 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
   double err_hi = ERR_HI[hi == 0.0];
 
   // Scaling factior = 2^(-xh_bits.get_exponent())
-  uint64_t s_u = (static_cast<uint64_t>(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];
+  int64_t s_u = static_cast<int64_t>(x_u & FPBits_t::EXP_MASK) -
+                (static_cast<int64_t>(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 = static_cast<uint64_t>(static_cast<int64_t>(x_u) - s_u);
+  uint64_t m_lo = (x_dd.lo != 0.0) ? FPBits_t(x_dd.lo).uintval() :
+                  static_cast<uint64_t>(cpp::bit_cast<int64_t>(x_dd.lo) - s_u);
+  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 eba65f56df7396..0d7b5ac80b2ff4 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 <stdint.h>
-
 using LlvmLibcLog1pTest = LIBC_NAMESPACE::testing::FPTest<double>;
 
 TEST_F(LlvmLibcLog1pTest, SpecialNumbers) {
@@ -26,6 +24,8 @@ 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 +36,21 @@ 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

[github-actions]

✅ With the latest revision this PR passed the C/C++ code formatter.

[michaelrj-google]

LGTM

[michaelrj-google]

since you flipped the subtraction in the calculation of s_u is this comment still accurate?

[lntue]

yes, it's still the same scaling factor, but we do subtract exponents (ldexp style) instead of multiply by 2^(-exp). I updated the comment to make it a bit clearer.