[libc-commits] [libc] [libc][math] Fix double-rounding errors with denormal inputs for hypot function. (PR #208247)

[via libc-commits]

llvmorg-github-actions[bot] wrote:


<!--LLVM PR SUMMARY COMMENT-->

@llvm/pr-subscribers-libc

Author: lntue

<details>
<summary>Changes</summary>

The issue was reported by Paul Zimmermann.

Solution: for denormal inputs, we use 1.0 + x trick after proper scaling to mimick rounding at denormal scale.

---
Full diff: https://github.com/llvm/llvm-project/pull/208247.diff


2 Files Affected:

- (modified) libc/src/__support/math/hypot.h (+71-4) 
- (modified) libc/test/src/math/smoke/hypot_test.cpp (+13) 


``````````diff
diff --git a/libc/src/__support/math/hypot.h b/libc/src/__support/math/hypot.h
index c92a700de1533..890f978e62aa5 100644
--- a/libc/src/__support/math/hypot.h
+++ b/libc/src/__support/math/hypot.h
@@ -22,6 +22,66 @@
 namespace LIBC_NAMESPACE_DECL {
 namespace math {
 
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+namespace hypot_internal {
+
+// Both input are denormals and non-zero.  We scale them up and down by the
+// inverse of the smallest normal number: 2^-1022, so that the correct "hidden"
+// bit position is now 1.  And by adding and subtracting 1 from the results, we
+// can emulate the rounding behavior in the denormal range.
+LIBC_INLINE double hypot_denorm(double a, double b) {
+  using fputil::DoubleDouble;
+  constexpr double SCALE = 0x1.0p1022;
+  constexpr double SCALE_BACK = 0x1.0p-1022;
+
+  a *= SCALE;
+  b *= SCALE;
+
+  // See the comments in the main function for the detail explanation of the
+  // computations.
+
+  // sum.hi + sum.lo ~ a^2 + b^2.
+  DoubleDouble a_sq = fputil::exact_mult(a, a);
+  DoubleDouble b_sq = fputil::exact_mult(b, b);
+  DoubleDouble sum = fputil::exact_add(a_sq.hi, b_sq.hi);
+  sum.lo += a_sq.lo + b_sq.lo;
+
+  // |sqrt(sum.hi) - r_hi| < 2^-52.
+  double r_hi = fputil::sqrt<double>(sum.hi);
+  // r_inv ~ 1 / (2 * r_hi)
+  double r_inv = 0.5 / r_hi;
+  // Adjust correction if needed.
+  DoubleDouble r_h{0.0, r_hi};
+  double correction = 0.0;
+  // 1 + rhi
+  if (r_hi < 1.0) {
+    r_h = fputil::exact_add(1.0, r_hi);
+    correction = 1.0;
+  }
+  // r_hi^2
+  DoubleDouble r_sq = fputil::exact_mult(r_hi, r_hi);
+  // (hi + lo - r_hi^2)
+  double num_lo = (sum.lo - r_sq.lo) - (r_sq.hi - sum.hi);
+  // (hi + lo - r_hi^2) / (2 * r_hi)
+  double r_lo = fputil::multiply_add(num_lo, r_inv, r_h.lo);
+
+  constexpr double ERR = 0x1.0p-100;
+
+  // Ziv's rounding test.
+  double upper = r_h.hi + (r_lo + ERR);
+  double lower = r_h.hi + (r_lo - ERR);
+
+  if (LIBC_LIKELY(upper == lower))
+    return (upper - correction) * SCALE_BACK;
+
+  return fputil::hypot(a * SCALE_BACK, b * SCALE_BACK);
+}
+
+} // namespace hypot_internal
+
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
 LIBC_INLINE double hypot(double x, double y) {
   using FPBits = fputil::FPBits<double>;
   using DoubleDouble = fputil::DoubleDouble;
@@ -80,10 +140,17 @@ LIBC_INLINE double hypot(double x, double y) {
   } else if (LIBC_UNLIKELY(b_e <= ((FPBits::EXP_BIAS - 500) << (32 - 11)))) {
     // The smaller magnitude is below 2^-500 (or 0), need to scale up to prevent
     // underflow when squaring.
-    if ((x == 0.0) || (y == 0.0)) {
-      double x_abs = FPBits(x_u & FPBits::EXP_SIG_MASK).get_val();
-      double y_abs = FPBits(y_u & FPBits::EXP_SIG_MASK).get_val();
-      return x_abs + y_abs;
+    if (LIBC_UNLIKELY(a_e < (1U << (32 - 11)))) {
+      // Larger input is denormal, extra care is needed to perform the Ziv's
+      // accuracy test correctly as double-rounding errors might happen.
+      if ((x == 0.0) || (y == 0.0)) {
+        double x_abs = FPBits(x_u & FPBits::EXP_SIG_MASK).get_val();
+        double y_abs = FPBits(y_u & FPBits::EXP_SIG_MASK).get_val();
+        return x_abs + y_abs;
+      }
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+      return hypot_internal::hypot_denorm(a, b);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
     }
     //  Any scaling factor > 2^((1072 + 52)/2) = 2^562 would work.
     scale = 0x1.0p600;
diff --git a/libc/test/src/math/smoke/hypot_test.cpp b/libc/test/src/math/smoke/hypot_test.cpp
index dd127909fafc5..b85e5ad456513 100644
--- a/libc/test/src/math/smoke/hypot_test.cpp
+++ b/libc/test/src/math/smoke/hypot_test.cpp
@@ -14,4 +14,17 @@ using LlvmLibcHypotTest = HypotTestTemplate<double>;
 
 TEST_F(LlvmLibcHypotTest, SpecialNumbers) {
   test_special_numbers(&LIBC_NAMESPACE::hypot);
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+  // Test denormal inputs.
+  EXPECT_FP_EQ(
+      0x0.c0bf7399534e3p-1022,
+      LIBC_NAMESPACE::hypot(0x0.2c2671b3c16b3p-1022, 0x0.bb9f8fecba9adp-1022));
+  EXPECT_FP_EQ(
+      0x1.6a09e667f3bcbp-1022,
+      LIBC_NAMESPACE::hypot(0x0.fffffffffffffp-1022, 0x0.fffffffffffffp-1022));
+  EXPECT_FP_EQ(
+      0x0.bffffb1b06483p-1022,
+      LIBC_NAMESPACE::hypot(0x0.603e52daf0bfdp-1022, -0x0.a622d0a9a433bp-1022));
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
 }

``````````

</details>


https://github.com/llvm/llvm-project/pull/208247