[libc][math] Add tests and fix some issues with FTZ/DAZ modes. (llvm#…

…113744)

libc/src/math/generic/atan2.cpp

@@ -230,8 +230,8 @@ LLVM_LIBC_FUNCTION(double, atan2, (double y, double x)) {
   if (LIBC_UNLIKELY(max_exp > 0x7ffU - 128U || min_exp < 128U)) {
     if (x_bits.is_nan() || y_bits.is_nan())
       return FPBits::quiet_nan().get_val();
-    unsigned x_except = x_abs == 0 ? 0 : (FPBits(x_abs).is_inf() ? 2 : 1);
-    unsigned y_except = y_abs == 0 ? 0 : (FPBits(y_abs).is_inf() ? 2 : 1);
+    unsigned x_except = x == 0.0 ? 0 : (FPBits(x_abs).is_inf() ? 2 : 1);
+    unsigned y_except = y == 0.0 ? 0 : (FPBits(y_abs).is_inf() ? 2 : 1);
 
     // Exceptional cases:
     //   EXCEPT[y_except][x_except][x_is_neg]

libc/src/math/generic/cbrt.cpp

@@ -151,9 +151,10 @@ LLVM_LIBC_FUNCTION(double, cbrt, (double x)) {
 
   if (LIBC_UNLIKELY(x_abs < FPBits::min_normal().uintval() ||
                     x_abs >= FPBits::inf().uintval())) {
-    if (x_abs == 0 || x_abs >= FPBits::inf().uintval())
+    if (x == 0.0 || x_abs >= FPBits::inf().uintval())
       // x is 0, Inf, or NaN.
-      return x;
+      // Make sure it works for FTZ/DAZ modes.
+      return static_cast<double>(x + x);
 
     // x is non-zero denormal number.
     // Normalize x.

libc/src/math/generic/cbrtf.cpp

@@ -93,9 +93,10 @@ LLVM_LIBC_FUNCTION(float, cbrtf, (float x)) {
   uint32_t x_abs = x_bits.uintval() & 0x7fff'ffff;
   uint32_t sign_bit = (x_bits.uintval() >> 31) << DoubleBits::EXP_LEN;
 
-  if (LIBC_UNLIKELY(x_abs == 0 || x_abs >= 0x7f80'0000)) {
+  if (LIBC_UNLIKELY(x == 0.0f || x_abs >= 0x7f80'0000)) {
     // x is 0, Inf, or NaN.
-    return x;
+    // Make sure it works for FTZ/DAZ modes.
+    return x + x;
   }
 
   double xd = static_cast<double>(x);

libc/src/math/generic/log.cpp

@@ -749,7 +749,7 @@ LLVM_LIBC_FUNCTION(double, log, (double x)) {
 
   if (LIBC_UNLIKELY(xbits.uintval() < FPBits_t::min_normal().uintval() ||
                     xbits.uintval() > FPBits_t::max_normal().uintval())) {
-    if (xbits.is_zero()) {
+    if (x == 0.0) {
       // return -Inf and raise FE_DIVBYZERO.
       fputil::set_errno_if_required(ERANGE);
       fputil::raise_except_if_required(FE_DIVBYZERO);

libc/src/math/generic/log10.cpp

@@ -751,7 +751,7 @@ LLVM_LIBC_FUNCTION(double, log10, (double x)) {
 
   if (LIBC_UNLIKELY(xbits.uintval() < FPBits_t::min_normal().uintval() ||
                     xbits.uintval() > FPBits_t::max_normal().uintval())) {
-    if (xbits.is_zero()) {
+    if (x == 0.0) {
       // return -Inf and raise FE_DIVBYZERO.
       fputil::set_errno_if_required(ERANGE);
       fputil::raise_except_if_required(FE_DIVBYZERO);

libc/src/math/generic/log10f.cpp

@@ -164,7 +164,7 @@ LLVM_LIBC_FUNCTION(float, log10f, (float x)) {
 
   if (LIBC_UNLIKELY(x_u < FPBits::min_normal().uintval() ||
                     x_u > FPBits::max_normal().uintval())) {
-    if (xbits.is_zero()) {
+    if (x == 0.0f) {
       // Return -inf and raise FE_DIVBYZERO
       fputil::set_errno_if_required(ERANGE);
       fputil::raise_except_if_required(FE_DIVBYZERO);

libc/src/math/generic/log1p.cpp

@@ -927,8 +927,8 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
       //   log(1 + x) = nextafter(x, -inf) for FE_DOWNWARD, or
       //                                       FE_TOWARDZERO and x > 0,
       //              = x                  otherwise.
-      if (LIBC_UNLIKELY(xbits.is_zero()))
-        return x;
+      if (x == 0.0)
+        return x + x; // Handle FTZ/DAZ correctly.
 
       volatile float tp = 1.0f;
       volatile float tn = -1.0f;
@@ -943,7 +943,7 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
         return FPBits_t(x_u + 1).get_val();
       }
 
-      return x;
+      return (x + x == 0.0) ? x + x : x;
     }
     x_dd = fputil::exact_add(1.0, x);
   }

libc/src/math/generic/log2.cpp

@@ -871,7 +871,7 @@ LLVM_LIBC_FUNCTION(double, log2, (double x)) {
 
   if (LIBC_UNLIKELY(xbits.uintval() < FPBits_t::min_normal().uintval() ||
                     xbits.uintval() > FPBits_t::max_normal().uintval())) {
-    if (xbits.is_zero()) {
+    if (x == 0.0) {
       // return -Inf and raise FE_DIVBYZERO.
       fputil::set_errno_if_required(ERANGE);
       fputil::raise_except_if_required(FE_DIVBYZERO);

libc/src/math/generic/log2f.cpp

@@ -72,7 +72,7 @@ LLVM_LIBC_FUNCTION(float, log2f, (float x)) {
   // Exceptional inputs.
   if (LIBC_UNLIKELY(x_u < FPBits::min_normal().uintval() ||
                     x_u > FPBits::max_normal().uintval())) {
-    if (xbits.is_zero()) {
+    if (x == 0.0f) {
       fputil::set_errno_if_required(ERANGE);
       fputil::raise_except_if_required(FE_DIVBYZERO);
       return FPBits::inf(Sign::NEG).get_val();

libc/src/math/generic/logf.cpp

@@ -82,7 +82,7 @@ LLVM_LIBC_FUNCTION(float, logf, (float x)) {
     }
     // Subnormal inputs.
     if (LIBC_UNLIKELY(x_u < FPBits::min_normal().uintval())) {
-      if (x_u == 0) {
+      if (x == 0.0f) {
         // Return -inf and raise FE_DIVBYZERO
         fputil::set_errno_if_required(ERANGE);
         fputil::raise_except_if_required(FE_DIVBYZERO);

libc/src/math/generic/pow.cpp

@@ -228,16 +228,18 @@ LLVM_LIBC_FUNCTION(double, pow, (double x, double y)) {
                     x_u >= FPBits::inf().uintval() ||
                     x_u < FPBits::min_normal().uintval())) {
     // Exceptional exponents.
-    switch (y_a) {
-    case 0: // y = +-0.0
+    if (y == 0.0)
       return 1.0;
+
+    switch (y_a) {
     case 0x3fe0'0000'0000'0000: { // y = +-0.5
       // TODO: speed up x^(-1/2) with rsqrt(x) when available.
-      if (LIBC_UNLIKELY(!y_sign && (x_u == FPBits::zero(Sign::NEG).uintval() ||
-                                    x_u == FPBits::inf(Sign::NEG).uintval()))) {
+      if (LIBC_UNLIKELY(
+              (x == 0.0 || x_u == FPBits::inf(Sign::NEG).uintval()))) {
         // pow(-0, 1/2) = +0
         // pow(-inf, 1/2) = +inf
-        return FPBits(x_abs).get_val();
+        // Make sure it works correctly for FTZ/DAZ.
+        return y_sign ? 1.0 / (x * x) : (x * x);
       }
       return y_sign ? (1.0 / fputil::sqrt<double>(x)) : fputil::sqrt<double>(x);
     }
@@ -269,7 +271,7 @@ LLVM_LIBC_FUNCTION(double, pow, (double x, double y)) {
           return 1.0;
         }
 
-        if (x_a == 0 && y_sign) {
+        if (x == 0.0 && y_sign) {
           // pow(+-0, -Inf) = +inf and raise FE_DIVBYZERO
           fputil::set_errno_if_required(EDOM);
           fputil::raise_except_if_required(FE_DIVBYZERO);
@@ -298,7 +300,7 @@ LLVM_LIBC_FUNCTION(double, pow, (double x, double y)) {
 
     // TODO: Speed things up with pow(2, y) = exp2(y) and pow(10, y) = exp10(y).
 
-    if (x_a == 0) {
+    if (x == 0.0) {
       bool out_is_neg = x_sign && is_odd_integer(y);
       if (y_sign) {
         // pow(0, negative number) = inf

libc/src/math/generic/powf.cpp

@@ -529,10 +529,10 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
   // Hence x^y will either overflow or underflow if x is not zero.
   if (LIBC_UNLIKELY((y_abs & 0x0007'ffff) == 0) || (y_abs > 0x4f170000)) {
     // Exceptional exponents.
-    switch (y_abs) {
-    case 0x0000'0000: { // y = +-0.0f
+    if (y == 0.0f)
       return 1.0f;
-    }
+
+    switch (y_abs) {
     case 0x7f80'0000: { // y = +-Inf
       if (x_abs > 0x7f80'0000) {
         // pow(NaN, +-Inf) = NaN
@@ -542,7 +542,7 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
         // pow(+-1, +-Inf) = 1.0f
         return 1.0f;
       }
-      if (x_abs == 0 && y_u == 0xff80'0000) {
+      if (x == 0.0f && y_u == 0xff80'0000) {
         // pow(+-0, -Inf) = +inf and raise FE_DIVBYZERO
         fputil::set_errno_if_required(EDOM);
         fputil::raise_except_if_required(FE_DIVBYZERO);
@@ -561,12 +561,15 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
       switch (y_u) {
       case 0x3f00'0000: // y = 0.5f
         // pow(x, 1/2) = sqrt(x)
-        if (LIBC_UNLIKELY(x_u == 0x8000'0000 || x_u == 0xff80'0000)) {
+        if (LIBC_UNLIKELY(x == 0.0f || x_u == 0xff80'0000)) {
           // pow(-0, 1/2) = +0
           // pow(-inf, 1/2) = +inf
-          return FloatBits(x_abs).get_val();
+          // Make sure it is correct for FTZ/DAZ.
+          return x * x;
         }
-        return fputil::sqrt<float>(x);
+        float r;
+        r = fputil::sqrt<float>(x);
+        return (FloatBits(r).uintval() != 0x8000'0000) ? r : 0.0f;
       case 0x3f80'0000: // y = 1.0f
         return x;
       case 0x4000'0000: // y = 2.0f
@@ -634,8 +637,7 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
 
     const bool x_is_neg = x_u >= FloatBits::SIGN_MASK;
 
-    switch (x_abs) {
-    case 0x0000'0000: { // x = +-0.0f
+    if (x == 0.0f) {
       const bool out_is_neg =
           x_is_neg && is_odd_integer(FloatBits(y_u).get_val());
       if (y_u > 0x8000'0000U) {
@@ -647,15 +649,16 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
       // pow(0, positive number) = 0
       return out_is_neg ? -0.0f : 0.0f;
     }
-    case 0x7f80'0000: { // x = +-Inf
+
+    if (x_abs == 0x7f80'0000) {
+      // x = +-Inf
       const bool out_is_neg =
           x_is_neg && is_odd_integer(FloatBits(y_u).get_val());
       if (y_u >= FloatBits::SIGN_MASK) {
         return out_is_neg ? -0.0f : 0.0f;
       }
       return FloatBits::inf(out_is_neg ? Sign::NEG : Sign::POS).get_val();
     }
-    }
 
     if (x_abs > 0x7f80'0000) {
       // x is NaN.

libc/src/math/generic/sin.cpp

@@ -50,7 +50,7 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
       if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 26)) {
         // Signed zeros.
         if (LIBC_UNLIKELY(x == 0.0))
-          return x;
+          return x + x; // Make sure it works with FTZ/DAZ.
 
 #ifdef LIBC_TARGET_CPU_HAS_FMA
         return fputil::multiply_add(x, -0x1.0p-54, x);

libc/src/math/generic/tan.cpp

@@ -138,7 +138,7 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) {
       if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) {
         // Signed zeros.
         if (LIBC_UNLIKELY(x == 0.0))
-          return x;
+          return x + x; // Make sure it works with FTZ/DAZ.
 
 #ifdef LIBC_TARGET_CPU_HAS_FMA
         return fputil::multiply_add(x, 0x1.0p-54, x);

libc/test/src/math/smoke/HypotTest.h

@@ -14,13 +14,11 @@
 #include "test/UnitTest/Test.h"
 
 template <typename T>
-class HypotTestTemplate : public LIBC_NAMESPACE::testing::Test {
-private:
+struct HypotTestTemplate : public LIBC_NAMESPACE::testing::Test {
   using Func = T (*)(T, T);
 
   DECLARE_SPECIAL_CONSTANTS(T)
 
-public:
   void test_special_numbers(Func func) {
     constexpr int N = 4;
     // Pythagorean triples.

libc/test/src/math/smoke/acosf_test.cpp

@@ -38,3 +38,27 @@ TEST_F(LlvmLibcAcosfTest, SpecialNumbers) {
   EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acosf(-2.0f));
   EXPECT_MATH_ERRNO(EDOM);
 }
+
+#ifdef LIBC_TEST_FTZ_DAZ
+
+using namespace LIBC_NAMESPACE::testing;
+
+TEST_F(LlvmLibcAcosfTest, FTZMode) {
+  ModifyMXCSR mxcsr(FTZ);
+
+  EXPECT_FP_EQ(0x1.921fb6p0f, LIBC_NAMESPACE::acosf(min_denormal));
+}
+
+TEST_F(LlvmLibcAcosfTest, DAZMode) {
+  ModifyMXCSR mxcsr(DAZ);
+
+  EXPECT_FP_EQ(0x1.921fb6p0f, LIBC_NAMESPACE::acosf(min_denormal));
+}
+
+TEST_F(LlvmLibcAcosfTest, FTZDAZMode) {
+  ModifyMXCSR mxcsr(FTZ | DAZ);
+
+  EXPECT_FP_EQ(0x1.921fb6p0f, LIBC_NAMESPACE::acosf(min_denormal));
+}
+
+#endif

libc/test/src/math/smoke/acoshf_test.cpp

@@ -35,3 +35,27 @@ TEST_F(LlvmLibcAcoshfTest, SpecialNumbers) {
   EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf(neg_inf));
   EXPECT_MATH_ERRNO(EDOM);
 }
+
+#ifdef LIBC_TEST_FTZ_DAZ
+
+using namespace LIBC_NAMESPACE::testing;
+
+TEST_F(LlvmLibcAcoshfTest, FTZMode) {
+  ModifyMXCSR mxcsr(FTZ);
+
+  EXPECT_FP_IS_NAN(LIBC_NAMESPACE::acoshf(min_denormal));
+}
+
+TEST_F(LlvmLibcAcoshfTest, DAZMode) {
+  ModifyMXCSR mxcsr(DAZ);
+
+  EXPECT_FP_IS_NAN(LIBC_NAMESPACE::acoshf(min_denormal));
+}
+
+TEST_F(LlvmLibcAcoshfTest, FTZDAZMode) {
+  ModifyMXCSR mxcsr(FTZ | DAZ);
+
+  EXPECT_FP_IS_NAN(LIBC_NAMESPACE::acoshf(min_denormal));
+}
+
+#endif

libc/test/src/math/smoke/asinf_test.cpp

@@ -41,3 +41,27 @@ TEST_F(LlvmLibcAsinfTest, SpecialNumbers) {
   EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::asinf(-2.0f));
   EXPECT_MATH_ERRNO(EDOM);
 }
+
+#ifdef LIBC_TEST_FTZ_DAZ
+
+using namespace LIBC_NAMESPACE::testing;
+
+TEST_F(LlvmLibcAsinfTest, FTZMode) {
+  ModifyMXCSR mxcsr(FTZ);
+
+  EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::asinf(min_denormal));
+}
+
+TEST_F(LlvmLibcAsinfTest, DAZMode) {
+  ModifyMXCSR mxcsr(DAZ);
+
+  EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::asinf(min_denormal));
+}
+
+TEST_F(LlvmLibcAsinfTest, FTZDAZMode) {
+  ModifyMXCSR mxcsr(FTZ | DAZ);
+
+  EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::asinf(min_denormal));
+}
+
+#endif

libc/test/src/math/smoke/asinhf_test.cpp

@@ -35,3 +35,27 @@ TEST_F(LlvmLibcAsinhfTest, SpecialNumbers) {
   EXPECT_FP_EQ_ALL_ROUNDING(neg_inf, LIBC_NAMESPACE::asinhf(neg_inf));
   EXPECT_MATH_ERRNO(0);
 }
+
+#ifdef LIBC_TEST_FTZ_DAZ
+
+using namespace LIBC_NAMESPACE::testing;
+
+TEST_F(LlvmLibcAsinhfTest, FTZMode) {
+  ModifyMXCSR mxcsr(FTZ);
+
+  EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::asinhf(min_denormal));
+}
+
+TEST_F(LlvmLibcAsinhfTest, DAZMode) {
+  ModifyMXCSR mxcsr(DAZ);
+
+  EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::asinhf(min_denormal));
+}
+
+TEST_F(LlvmLibcAsinhfTest, FTZDAZMode) {
+  ModifyMXCSR mxcsr(FTZ | DAZ);
+
+  EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::asinhf(min_denormal));
+}
+
+#endif