[libc] Fix bugs with negative and mixed normal/denormal inputs in hypot implementation.

[lntue]

Summary

Fix bugs with negative and mixed normal/denormal inputs in hypot implementation.

[lntue]

[libc] Fix bugs with negative and mixed normal/denormal inputs in hypot implementation.

[michaelrj]

libc/utils/MPFRWrapper/MPFRUtils.cpp

315–318 ↗

(On Diff #390746)

based on my testing, the tests still pass if you don't update this file, could you add something to explain these changes, and possibly split them off into their own patch?

[lntue]

libc/utils/MPFRWrapper/MPFRUtils.cpp

315–318 ↗

(On Diff #390746)

This is because without this change, the ULPs are computed off by a factor of 2^(-mantissaWidth), and hence all the tests still pass. As in the comments in line 287-294, the factor that we divide is eps(...) which should be 2^(exponent - mantissaWidth) for normal numbers, and 2^(1 - exponentBias - mantissaWidth) for denormal numbers.

With changes in this file, updated HypotTests without the fix in Hypot.h shall fail. In any case, I split the change in the ULP function to https://reviews.llvm.org/D114726

[lntue]

[libc] Fix bugs with negative and mixed normal/denormal inputs in hypot implementation.

[lntue]

[libc] Fix bugs with negative and mixed normal/denormal inputs in hypot implementation.

[lntue]

[libc] Fix bugs with negative and mixed normal/denormal inputs in hypot implementation.

[lntue]

[libc] Fix bugs with negative and mixed normal/denormal inputs in hypot implementation.

libc/src/__support/FPUtil/Hypot.h

if ((x_bits.getUnbiasedExponent() >=
     return abs(x);
   } else if ((y_bits.getUnbiasedExponent() >=
               x_bits.getUnbiasedExponent() + MantissaWidth<T>::value + 2) ||
              (x == 0)) {
     y_bits.setSign(0);
     return abs(y);
   }
 
-  if (x >= y) {
+  if (abs(x) >= abs(y)) {
     a_exp = x_bits.getUnbiasedExponent();
     a_mant = x_bits.getMantissa();
     b_exp = y_bits.getUnbiasedExponent();
     b_mant = y_bits.getMantissa();
   } else {
     a_exp = y_bits.getUnbiasedExponent();
     a_mant = y_bits.getMantissa();
     b_exp = x_bits.getUnbiasedExponent();
static inline T hypot(T x, T y) {
   UIntType leading_one;
   int y_mant_width;
   if (a_exp != 0) {
     leading_one = one;
     a_mant |= one;
     y_mant_width = MantissaWidth<T>::value + 1;
   } else {
     leading_one = internal::findLeadingOne(a_mant, y_mant_width);
+    a_exp = 1;
   }
 
   if (b_exp != 0) {
     b_mant |= one;
+  } else {
+    b_exp = 1;
   }
 
   a_mant_sq = static_cast<DUIntType>(a_mant) * a_mant;
   b_mant_sq = static_cast<DUIntType>(b_mant) * b_mant;
 
   // At this point, a_exp >= b_exp > a_exp - 25, so in order to line up aSqMant
   // and bSqMant, we need to shift bSqMant to the right by (a_exp - b_exp) bits.
   // But before that, remember to store the losing bits to sticky.

libc/test/src/math/HypotTest.h

void testSpecialNumbers(Func func) {
     EXPECT_FP_EQ(func(nan, zero), nan);
     EXPECT_FP_EQ(func(negZero, nan), nan);
 
     EXPECT_FP_EQ(func(negZero, zero), zero);
   }
 
   void testSubnormalRange(Func func) {
     constexpr UIntType count = 1000001;
-    constexpr UIntType step =
-        (FPBits::maxSubnormal - FPBits::minSubnormal) / count;
-    for (UIntType v = FPBits::minSubnormal, w = FPBits::maxSubnormal;
-         v <= FPBits::maxSubnormal && w >= FPBits::minSubnormal;
-         v += step, w -= step) {
+    for (unsigned scale = 0; scale < 4; ++scale) {
+      UIntType maxValue = FPBits::maxSubnormal << scale;
+      UIntType step = (maxValue - FPBits::minSubnormal) / count;
+      for (int signs = 0; signs < 4; ++signs) {
+        for (UIntType v = FPBits::minSubnormal, w = maxValue;
+             v <= maxValue && w >= FPBits::minSubnormal; v += step, w -= step) {
           T x = T(FPBits(v)), y = T(FPBits(w));
+          if (signs % 2 == 1) {
+            x = -x;
+          }
+          if (signs >= 2) {
+            y = -y;
+          }
+
           T result = func(x, y);
           mpfr::BinaryInput<T> input{x, y};
           ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5);
         }
       }
+    }
+  }
 
   void testNormalRange(Func func) {
     constexpr UIntType count = 1000001;
     constexpr UIntType step = (FPBits::maxNormal - FPBits::minNormal) / count;
+    for (int signs = 0; signs < 4; ++signs) {
       for (UIntType v = FPBits::minNormal, w = FPBits::maxNormal;
            v <= FPBits::maxNormal && w >= FPBits::minNormal;
            v += step, w -= step) {
         T x = T(FPBits(v)), y = T(FPBits(w));
+        if (signs % 2 == 1) {
+          x = -x;
+        }
+        if (signs >= 2) {
+          y = -y;
+        }
+
         T result = func(x, y);
         mpfr::BinaryInput<T> input{x, y};
         ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5);
       }
     }
+  }
 };
 
 #endif // LLVM_LIBC_TEST_SRC_MATH_HYPOTTEST_H