Fix when m_lo is underflowed.

Author: lntue

libc/src/math/generic/log1p.cpp

@@ -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<int>(((x_u & FRACTION_MASK) + (1ULL << (FRACTION_LEN - 8))) >>
-                       (FRACTION_LEN - 7));
+  int idx = static_cast<int>((xhi_frac + (1ULL << (FRACTION_LEN - 8))) >>
+                             (FRACTION_LEN - 7));
   int x_e = xhi_bits.get_exponent() + (idx >> 7);
   double e_x = static_cast<double>(x_e);
 
@@ -981,11 +980,13 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
   //   1 <= m_dd.hi < 2
   //   |m_dd.lo| < 2^-52.
   // This is exact.
-  uint64_t m_hi = static_cast<uint64_t>(static_cast<int64_t>(x_u) - s_u);
+  uint64_t m_hi = FPBits_t::one().uintval() | xhi_frac;
+
   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);
+      FPBits_t(x_dd.lo).abs().get_val() > x_dd.hi * 0x1.0p-127
+          ? static_cast<uint64_t>(cpp::bit_cast<int64_t>(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: