[libc-commits] [libc] [libc][math] Fix log1p SEGV with large inputs when FTZ/DAZ flags are set. (PR #115541)
via libc-commits
libc-commits at lists.llvm.org
Fri Nov 8 19:30:44 PST 2024
https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/115541
>From a3b6177c3f2197ac5cf24d566ea77ca56476d3de Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 8 Nov 2024 20:20:45 +0000
Subject: [PATCH 1/5] [libc][math] Fix log1p SEGV with large inputs when
FTZ/DAZ flags are set.
---
libc/src/math/generic/log1p.cpp | 14 +++++++-------
libc/test/src/math/smoke/log1p_test.cpp | 7 +++++--
2 files changed, 12 insertions(+), 9 deletions(-)
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
>From e34408dc36d2c96c91658d128455faa8baa73f03 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 8 Nov 2024 20:38:06 +0000
Subject: [PATCH 2/5] format.
---
libc/src/math/generic/log1p.cpp | 8 +++++---
libc/test/src/math/smoke/log1p_test.cpp | 12 ++++++++----
2 files changed, 13 insertions(+), 7 deletions(-)
diff --git a/libc/src/math/generic/log1p.cpp b/libc/src/math/generic/log1p.cpp
index ec1f2bc00fd81f..38f3bed05d5796 100644
--- a/libc/src/math/generic/log1p.cpp
+++ b/libc/src/math/generic/log1p.cpp
@@ -823,7 +823,7 @@ constexpr Float128 BIG_COEFFS[4]{
};
[[maybe_unused]] LIBC_INLINE double log1p_accurate(int e_x, int index,
- fputil::DoubleDouble m_x) {
+ fputil::DoubleDouble m_x) {
Float128 e_x_f128(static_cast<float>(e_x));
Float128 sum = fputil::quick_mul(LOG_2, e_x_f128);
sum = fputil::quick_add(sum, LOG_R1[index]);
@@ -982,8 +982,10 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
// |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_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);
+ 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:
diff --git a/libc/test/src/math/smoke/log1p_test.cpp b/libc/test/src/math/smoke/log1p_test.cpp
index 0d7b5ac80b2ff4..b98c0f26a8bcae 100644
--- a/libc/test/src/math/smoke/log1p_test.cpp
+++ b/libc/test/src/math/smoke/log1p_test.cpp
@@ -25,7 +25,8 @@ TEST_F(LlvmLibcLog1pTest, SpecialNumbers) {
EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, LIBC_NAMESPACE::log1p(-1.0),
FE_DIVBYZERO);
- EXPECT_FP_EQ(0x1.62c829bf8fd9dp9, LIBC_NAMESPACE::log1p(0x1.9b536cac3a09dp1023));
+ EXPECT_FP_EQ(0x1.62c829bf8fd9dp9,
+ LIBC_NAMESPACE::log1p(0x1.9b536cac3a09dp1023));
}
#ifdef LIBC_TEST_FTZ_DAZ
@@ -36,21 +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));
+ 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));
+ 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));
+ EXPECT_FP_EQ(0x1.62c829bf8fd9dp9,
+ LIBC_NAMESPACE::log1p(0x1.9b536cac3a09dp1023));
}
#endif
>From ba3355689a6339c02c45cff8b9e4762d827b102b Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 8 Nov 2024 20:47:36 +0000
Subject: [PATCH 3/5] Fix typo.
---
libc/src/math/generic/log1p.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/src/math/generic/log1p.cpp b/libc/src/math/generic/log1p.cpp
index 38f3bed05d5796..31f51fa9b33bde 100644
--- a/libc/src/math/generic/log1p.cpp
+++ b/libc/src/math/generic/log1p.cpp
@@ -974,7 +974,7 @@ 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())
+ // Scaling factor = 2^(-xh_bits.get_exponent())
int64_t s_u = static_cast<int64_t>(x_u & FPBits_t::EXP_MASK) -
(static_cast<int64_t>(EXP_BIAS) << FRACTION_LEN);
// Normalize arguments:
>From ab057bdb57ba30a99a42b0d6ef8daa8fb4b60315 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 8 Nov 2024 20:49:01 +0000
Subject: [PATCH 4/5] Update comment.
---
libc/src/math/generic/log1p.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/src/math/generic/log1p.cpp b/libc/src/math/generic/log1p.cpp
index 31f51fa9b33bde..a7721edcf9e763 100644
--- a/libc/src/math/generic/log1p.cpp
+++ b/libc/src/math/generic/log1p.cpp
@@ -974,7 +974,7 @@ 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 factor = 2^(-xh_bits.get_exponent())
+ // Scale x_dd by 2^(-xh_bits.get_exponent()).
int64_t s_u = static_cast<int64_t>(x_u & FPBits_t::EXP_MASK) -
(static_cast<int64_t>(EXP_BIAS) << FRACTION_LEN);
// Normalize arguments:
>From 4bf657fa86a10fa5efb7fd6edffd47a37a1be101 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Sat, 9 Nov 2024 03:29:37 +0000
Subject: [PATCH 5/5] Fix when m_lo is underflowed.
---
libc/src/math/generic/log1p.cpp | 17 +++++++++--------
1 file changed, 9 insertions(+), 8 deletions(-)
diff --git a/libc/src/math/generic/log1p.cpp b/libc/src/math/generic/log1p.cpp
index a7721edcf9e763..b9c58b843a2409 100644
--- a/libc/src/math/generic/log1p.cpp
+++ b/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:
More information about the libc-commits
mailing list