Skip to content

Commit

Permalink
[libc][math] Fix signed zeros for erff. (llvm#97742)
Browse files Browse the repository at this point in the history
The inexact exception flag was raised for the exact cases of signed
zeros. This was reported by Paul Zimmermann using the CORE-MATH test
suites.
  • Loading branch information
lntue committed Jul 5, 2024
1 parent d4216b5 commit 9f2215a
Showing 1 changed file with 15 additions and 9 deletions.
24 changes: 15 additions & 9 deletions libc/src/math/generic/erff.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -127,15 +127,6 @@ LLVM_LIBC_FUNCTION(float, erff, (float x)) {
uint32_t x_u = xbits.uintval();
uint32_t x_abs = x_u & 0x7fff'ffffU;

// Exceptional values
if (LIBC_UNLIKELY(x_abs == 0x3f65'9229U)) // |x| = 0x1.cb2452p-1f
return x < 0.0f ? fputil::round_result_slightly_down(-0x1.972ea8p-1f)
: fputil::round_result_slightly_up(0x1.972ea8p-1f);
if (LIBC_UNLIKELY(x_abs == 0x4004'1e6aU)) // |x| = 0x1.083cd4p+1f
return x < 0.0f ? fputil::round_result_slightly_down(-0x1.fe3462p-1f)
: fputil::round_result_slightly_up(0x1.fe3462p-1f);

// if (LIBC_UNLIKELY(x_abs > 0x407a'd444U)) {
if (LIBC_UNLIKELY(x_abs >= 0x4080'0000U)) {
const float ONE[2] = {1.0f, -1.0f};
const float SMALL[2] = {-0x1.0p-25f, 0x1.0p-25f};
Expand All @@ -149,6 +140,21 @@ LLVM_LIBC_FUNCTION(float, erff, (float x)) {
return ONE[sign] + SMALL[sign];
}

// Exceptional mask = common 0 bits of 2 exceptional values.
constexpr uint32_t EXCEPT_MASK = 0x809a'6184U;

if (LIBC_UNLIKELY((x_abs & EXCEPT_MASK) == 0)) {
// Exceptional values
if (LIBC_UNLIKELY(x_abs == 0x3f65'9229U)) // |x| = 0x1.cb2452p-1f
return x < 0.0f ? fputil::round_result_slightly_down(-0x1.972ea8p-1f)
: fputil::round_result_slightly_up(0x1.972ea8p-1f);
if (LIBC_UNLIKELY(x_abs == 0x4004'1e6aU)) // |x| = 0x1.083cd4p+1f
return x < 0.0f ? fputil::round_result_slightly_down(-0x1.fe3462p-1f)
: fputil::round_result_slightly_up(0x1.fe3462p-1f);
if (x_abs == 0U)
return x;
}

// Polynomial approximation:
// erf(x) ~ x * (c0 + c1 * x^2 + c2 * x^4 + ... + c7 * x^14)
double xd = static_cast<double>(x);
Expand Down

0 comments on commit 9f2215a

Please sign in to comment.