[libc-commits] [libc] [libc][math] Fix double-rounding errors with denormal inputs for hypot function. (PR #208247)
via libc-commits
libc-commits at lists.llvm.org
Thu Jul 9 06:54:54 PDT 2026
https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/208247
>From 8b4453eaa18723dc06db7190dd97c27d57be38f6 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Wed, 8 Jul 2026 15:57:23 +0000
Subject: [PATCH 1/5] [libc][math] Fix double-rounding errors with denormal
inputs for hypot function.
---
libc/src/__support/math/hypot.h | 75 +++++++++++++++++++++++--
libc/test/src/math/smoke/hypot_test.cpp | 13 +++++
2 files changed, 84 insertions(+), 4 deletions(-)
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
}
>From ce26ec954e4c9efe3248a67feac993dd650a4572 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Wed, 8 Jul 2026 18:02:13 +0000
Subject: [PATCH 2/5] Update comment and error bound.
---
libc/src/__support/math/hypot.h | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/libc/src/__support/math/hypot.h b/libc/src/__support/math/hypot.h
index 890f978e62aa5..a919638a9e939 100644
--- a/libc/src/__support/math/hypot.h
+++ b/libc/src/__support/math/hypot.h
@@ -54,8 +54,9 @@ LIBC_INLINE double hypot_denorm(double a, double b) {
// Adjust correction if needed.
DoubleDouble r_h{0.0, r_hi};
double correction = 0.0;
- // 1 + rhi
if (r_hi < 1.0) {
+ // When r_hi < 1, the output is denormal. We mimick rounding in denormal
+ // range with 1.0 + r_hi.
r_h = fputil::exact_add(1.0, r_hi);
correction = 1.0;
}
@@ -66,7 +67,7 @@ LIBC_INLINE double hypot_denorm(double a, double b) {
// (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;
+ constexpr double ERR = 0x1.0p-102;
// Ziv's rounding test.
double upper = r_h.hi + (r_lo + ERR);
>From 8c437d91174f8e7e9c8b4bb8faa3098dbd119fc6 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Wed, 8 Jul 2026 22:29:15 +0000
Subject: [PATCH 3/5] Fix overflow and inexact exceptions.
---
libc/src/__support/FPUtil/Hypot.h | 6 ++++--
libc/src/__support/math/hypot.h | 17 +++++++++++++++++
2 files changed, 21 insertions(+), 2 deletions(-)
diff --git a/libc/src/__support/FPUtil/Hypot.h b/libc/src/__support/FPUtil/Hypot.h
index fdbd7b4acf088..bd20fdfbb719b 100644
--- a/libc/src/__support/FPUtil/Hypot.h
+++ b/libc/src/__support/FPUtil/Hypot.h
@@ -283,8 +283,10 @@ LIBC_INLINE T hypot(T x, T y) {
y_new |= static_cast<StorageType>(out_exp) << FPBits_t::FRACTION_LEN;
- if (!(round_bit || sticky_bits || (r != 0)))
- fputil::clear_except_if_required(FE_INEXACT);
+ // TODO: We should only clear FE_INEXACT except if it's not set at the
+ // start of the function.
+ // if (!(round_bit || sticky_bits || (r != 0)))
+ // fputil::clear_except_if_required(FE_INEXACT);
return cpp::bit_cast<T>(y_new);
}
diff --git a/libc/src/__support/math/hypot.h b/libc/src/__support/math/hypot.h
index a919638a9e939..7688013b195d1 100644
--- a/libc/src/__support/math/hypot.h
+++ b/libc/src/__support/math/hypot.h
@@ -138,6 +138,23 @@ LIBC_INLINE double hypot(double x, double y) {
scale_back = 0x1.0p600;
a *= scale;
b *= scale;
+ // Check for overflow to raise them correctly.
+#if !defined(LIBC_MATH_HAS_NO_EXCEPT)
+ // No overflow when calculating a^2 + b^2.
+ double asq = a * a;
+ double bsq = b * b;
+ double sumsq = asq + bsq;
+ // Overflow happens when:
+ // 2^600 * sqrt(a^2 + b^2) >= 2^1023 * (2 - 2^-53)
+ // Which is equivalent to:
+ // sqrt(a^2 + b^2) >= 2^424 * (1 - 2^-54).
+ // Square both sides:
+ // a^2 + b^2 >= 2^848 * (1 - 2^-53 + 2^-108).
+ // For a fast sufficient condition that can be done in double precision:
+ // a^2 + b^2 >= 2^848.
+ if (sumsq >= 0x1.0p848)
+ return sumsq * scale_back;
+#endif // !LIBC_MATH_HAS_NO_EXCEPT
} 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.
>From 80ab28a5ebe74bd99512df7e6c96a0ccddc2a264 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Thu, 9 Jul 2026 00:57:58 +0000
Subject: [PATCH 4/5] Raise overflow and underflow exceptions properly.
---
libc/src/__support/math/hypot.h | 28 ++++++++++++++++++++++++----
1 file changed, 24 insertions(+), 4 deletions(-)
diff --git a/libc/src/__support/math/hypot.h b/libc/src/__support/math/hypot.h
index 7688013b195d1..3b33503506cf6 100644
--- a/libc/src/__support/math/hypot.h
+++ b/libc/src/__support/math/hypot.h
@@ -73,8 +73,20 @@ LIBC_INLINE double hypot_denorm(double a, double b) {
double upper = r_h.hi + (r_lo + ERR);
double lower = r_h.hi + (r_lo - ERR);
- if (LIBC_LIKELY(upper == lower))
+ if (LIBC_LIKELY(upper == lower)) {
+#ifdef LIBC_MATH_HAS_NO_EXCEPT
return (upper - correction) * SCALE_BACK;
+#else
+ // Check to raise underflow correctly.
+ DoubleDouble r = fputil::exact_add(r_h.hi, r_lo);
+ r.hi -= correction;
+ // Raise underflow if needed:
+ if ((r.hi < 1.0 && r.lo != 0.0) || (r.hi == 1.0 && r.lo < 0.0))
+ fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
+
+ return r.hi * SCALE_BACK;
+#endif // LIBC_MATH_HAS_NO_EXCEPT
+ }
return fputil::hypot(a * SCALE_BACK, b * SCALE_BACK);
}
@@ -133,12 +145,20 @@ LIBC_INLINE double hypot(double x, double y) {
return x;
return y;
}
+
+ // Check the exponent gap here so that all the follow up pre-scaling and
+ // overflow check won't generate spurious underflow exceptions.
+ if (LIBC_UNLIKELY(a_e - b_e >= (54U << (32 - 11)))) {
+ 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;
+ }
// Any scaling factor < 2^(-1024/2) = 2^-512 would work.
scale = 0x1.0p-600;
scale_back = 0x1.0p600;
a *= scale;
b *= scale;
- // Check for overflow to raise them correctly.
+ // Check for overflow to raise the exception correctly.
#if !defined(LIBC_MATH_HAS_NO_EXCEPT)
// No overflow when calculating a^2 + b^2.
double asq = a * a;
@@ -155,8 +175,8 @@ LIBC_INLINE double hypot(double x, double y) {
if (sumsq >= 0x1.0p848)
return sumsq * scale_back;
#endif // !LIBC_MATH_HAS_NO_EXCEPT
- } 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
+ } else if (LIBC_UNLIKELY(b_e <= ((FPBits::EXP_BIAS - 400) << (32 - 11)))) {
+ // The smaller magnitude is below 2^-400 (or 0), need to scale up to prevent
// underflow when squaring.
if (LIBC_UNLIKELY(a_e < (1U << (32 - 11)))) {
// Larger input is denormal, extra care is needed to perform the Ziv's
>From 3127315d6006c41a067a5521069c3a42a57f24c0 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Thu, 9 Jul 2026 13:54:22 +0000
Subject: [PATCH 5/5] Change to raise underflow after rounding.
---
libc/src/__support/math/hypot.h | 4 +++-
1 file changed, 3 insertions(+), 1 deletion(-)
diff --git a/libc/src/__support/math/hypot.h b/libc/src/__support/math/hypot.h
index 3b33503506cf6..aa17d49e76b6d 100644
--- a/libc/src/__support/math/hypot.h
+++ b/libc/src/__support/math/hypot.h
@@ -19,6 +19,8 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
+// #include <iostream>
+
namespace LIBC_NAMESPACE_DECL {
namespace math {
@@ -81,7 +83,7 @@ LIBC_INLINE double hypot_denorm(double a, double b) {
DoubleDouble r = fputil::exact_add(r_h.hi, r_lo);
r.hi -= correction;
// Raise underflow if needed:
- if ((r.hi < 1.0 && r.lo != 0.0) || (r.hi == 1.0 && r.lo < 0.0))
+ if (r.hi < 1.0 && r.lo != 0.0)
fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
return r.hi * SCALE_BACK;
More information about the libc-commits
mailing list