From 9f2215ae5517a3d6e1d59e597bff83b821b96c7e Mon Sep 17 00:00:00 2001
From: lntue <35648136+lntue@users.noreply.github.com>
Date: Fri, 5 Jul 2024 09:54:37 -0400
Subject: [PATCH] [libc][math] Fix signed zeros for erff. (#97742)

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.
---
 libc/src/math/generic/erff.cpp | 24 +++++++++++++++---------
 1 file changed, 15 insertions(+), 9 deletions(-)

diff --git a/libc/src/math/generic/erff.cpp b/libc/src/math/generic/erff.cpp
index f120d5646e0439..aa7baffc7815e9 100644
--- a/libc/src/math/generic/erff.cpp
+++ b/libc/src/math/generic/erff.cpp
@@ -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};
@@ -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);