[libcxx-commits] [libcxx] [libc++][math] Fix undue overflowing of `std::hypot(x, y, z)` (PR #93350)
via libcxx-commits
libcxx-commits at lists.llvm.org
Fri May 24 14:28:24 PDT 2024
llvmbot wrote:
<!--LLVM PR SUMMARY COMMENT-->
@llvm/pr-subscribers-libcxx
Author: None (PaulXiCao)
<details>
<summary>Changes</summary>
Closes #<!-- -->92782.
The 3-dimentionsional `std::hypot(x,y,z)` was sub-optimally implemented. This lead to possible over-/underflows in (intermediate) results which can be circumvented by this proposed change.
The idea is to to scale the arguments (see linked issue for full discussion).
Tests have been added for problematic over- and underflows.
---
Full diff: https://github.com/llvm/llvm-project/pull/93350.diff
3 Files Affected:
- (modified) libcxx/include/__math/hypot.h (+78)
- (modified) libcxx/include/cmath (+1-24)
- (modified) libcxx/test/std/numerics/c.math/cmath.pass.cpp (+62-15)
``````````diff
diff --git a/libcxx/include/__math/hypot.h b/libcxx/include/__math/hypot.h
index 8e2c0f5b6b6d8..6b68e7bfed341 100644
--- a/libcxx/include/__math/hypot.h
+++ b/libcxx/include/__math/hypot.h
@@ -9,11 +9,16 @@
#ifndef _LIBCPP___MATH_HYPOT_H
#define _LIBCPP___MATH_HYPOT_H
+#include <__algorithm/max.h>
#include <__config>
+#include <__math/abs.h>
+#include <__math/roots.h>
#include <__type_traits/enable_if.h>
#include <__type_traits/is_arithmetic.h>
#include <__type_traits/is_same.h>
#include <__type_traits/promote.h>
+#include <array>
+#include <limits>
#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
# pragma GCC system_header
@@ -41,6 +46,79 @@ inline _LIBCPP_HIDE_FROM_ABI typename __promote<_A1, _A2>::type hypot(_A1 __x, _
return __math::hypot((__result_type)__x, (__result_type)__y);
}
+#if _LIBCPP_STD_VER >= 17
+template <class _Real>
+struct __hypot_factors {
+ _Real __threshold;
+ _Real __scale_xyz;
+ _Real __scale_M;
+};
+
+// returns [underflow_factors, overflow_factors]
+template <class _Real>
+std::array<__hypot_factors<_Real>, 2> __create_factors() {
+ static_assert(std::numeric_limits<_Real>::is_iec559);
+
+ __hypot_factors<_Real> __underflow, __overflow;
+ if constexpr (std::is_same_v<_Real, float>) {
+ static_assert(-125 == std::numeric_limits<_Real>::min_exponent);
+ static_assert(+128 == std::numeric_limits<_Real>::max_exponent);
+ __underflow = __hypot_factors<_Real>{0x1.0p-62f, 0x1.0p70f, 0x1.0p-70f};
+ __overflow = __hypot_factors<_Real>{0x1.0p62f, 0x1.0p-70f, 0x1.0p+70f};
+ } else if constexpr (std::is_same_v<_Real, double>) {
+ static_assert(-1021 == std::numeric_limits<_Real>::min_exponent);
+ static_assert(+1024 == std::numeric_limits<_Real>::max_exponent);
+ __underflow = __hypot_factors<_Real>{0x1.0p-510, 0x1.0p600, 0x1.0p-600};
+ __overflow = __hypot_factors<_Real>{0x1.0p510, 0x1.0p-600, 0x1.0p+600};
+ } else { // long double
+ static_assert(std::is_same_v<_Real, long double>);
+ static_assert(-16'381 == std::numeric_limits<_Real>::min_exponent);
+ static_assert(+16'384 == std::numeric_limits<_Real>::max_exponent);
+ __underflow = __hypot_factors<_Real>{0x1.0p-8'190l, 0x1.0p9'000l, 0x1.0p-9'000l};
+ __overflow = __hypot_factors<_Real>{0x1.0p8'190l, 0x1.0p-9'000l, 0x1.0p+9'000l};
+ }
+ return {__underflow, __overflow};
+}
+
+template <class _Real>
+_Real __hypot(_Real __x, _Real __y, _Real __z) {
+ const auto [__underflow, __overflow] = __create_factors<_Real>();
+ _Real __M = std::max({__math::fabs(__x), __math::fabs(__y), __math::fabs(__z)});
+ if (__M > __overflow.__threshold) { // x*x + y*y + z*z might overflow
+ __x *= __overflow.__scale_xyz;
+ __y *= __overflow.__scale_xyz;
+ __z *= __overflow.__scale_xyz;
+ __M = __overflow.__scale_M;
+ } else if (__M < __underflow.__threshold) { // x*x + y*y + z*z might underflow
+ __x *= __underflow.__scale_xyz;
+ __y *= __underflow.__scale_xyz;
+ __z *= __underflow.__scale_xyz;
+ __M = __underflow.__scale_M;
+ } else
+ __M = 1;
+ return __M * __math::sqrt(__x * __x + __y * __y + __z * __z);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI float hypot(float __x, float __y, float __z) { return __hypot(__x, __y, __z); }
+
+inline _LIBCPP_HIDE_FROM_ABI double hypot(double __x, double __y, double __z) { return __hypot(__x, __y, __z); }
+
+inline _LIBCPP_HIDE_FROM_ABI long double hypot(long double __x, long double __y, long double __z) {
+ return __hypot(__x, __y, __z);
+}
+
+template <class _A1,
+ class _A2,
+ class _A3,
+ std::enable_if_t< is_arithmetic_v<_A1> && is_arithmetic_v<_A2> && is_arithmetic_v<_A3>, int> = 0 >
+_LIBCPP_HIDE_FROM_ABI typename __promote<_A1, _A2, _A3>::type hypot(_A1 __x, _A2 __y, _A3 __z) _NOEXCEPT {
+ using __result_type = typename __promote<_A1, _A2, _A3>::type;
+ static_assert(!(
+ std::is_same_v<_A1, __result_type> && std::is_same_v<_A2, __result_type> && std::is_same_v<_A3, __result_type>));
+ return __hypot(static_cast<__result_type>(__x), static_cast<__result_type>(__y), static_cast<__result_type>(__z));
+}
+#endif
+
} // namespace __math
_LIBCPP_END_NAMESPACE_STD
diff --git a/libcxx/include/cmath b/libcxx/include/cmath
index dd194bbb55896..49f35caa7b74e 100644
--- a/libcxx/include/cmath
+++ b/libcxx/include/cmath
@@ -305,6 +305,7 @@ constexpr long double lerp(long double a, long double b, long double t) noexcept
*/
#include <__config>
+#include <__math/hypot.h>
#include <__type_traits/enable_if.h>
#include <__type_traits/is_arithmetic.h>
#include <__type_traits/is_constant_evaluated.h>
@@ -544,30 +545,6 @@ using ::scalbnl _LIBCPP_USING_IF_EXISTS;
using ::tgammal _LIBCPP_USING_IF_EXISTS;
using ::truncl _LIBCPP_USING_IF_EXISTS;
-#if _LIBCPP_STD_VER >= 17
-inline _LIBCPP_HIDE_FROM_ABI float hypot(float __x, float __y, float __z) {
- return sqrt(__x * __x + __y * __y + __z * __z);
-}
-inline _LIBCPP_HIDE_FROM_ABI double hypot(double __x, double __y, double __z) {
- return sqrt(__x * __x + __y * __y + __z * __z);
-}
-inline _LIBCPP_HIDE_FROM_ABI long double hypot(long double __x, long double __y, long double __z) {
- return sqrt(__x * __x + __y * __y + __z * __z);
-}
-
-template <class _A1, class _A2, class _A3>
-inline _LIBCPP_HIDE_FROM_ABI
- typename enable_if_t< is_arithmetic<_A1>::value && is_arithmetic<_A2>::value && is_arithmetic<_A3>::value,
- __promote<_A1, _A2, _A3> >::type
- hypot(_A1 __lcpp_x, _A2 __lcpp_y, _A3 __lcpp_z) _NOEXCEPT {
- typedef typename __promote<_A1, _A2, _A3>::type __result_type;
- static_assert((!(is_same<_A1, __result_type>::value && is_same<_A2, __result_type>::value &&
- is_same<_A3, __result_type>::value)),
- "");
- return std::hypot((__result_type)__lcpp_x, (__result_type)__lcpp_y, (__result_type)__lcpp_z);
-}
-#endif
-
template <class _A1, __enable_if_t<is_floating_point<_A1>::value, int> = 0>
_LIBCPP_HIDE_FROM_ABI _LIBCPP_CONSTEXPR bool __constexpr_isnan(_A1 __lcpp_x) _NOEXCEPT {
#if __has_builtin(__builtin_isnan)
diff --git a/libcxx/test/std/numerics/c.math/cmath.pass.cpp b/libcxx/test/std/numerics/c.math/cmath.pass.cpp
index 9379084499792..544c197ceb1af 100644
--- a/libcxx/test/std/numerics/c.math/cmath.pass.cpp
+++ b/libcxx/test/std/numerics/c.math/cmath.pass.cpp
@@ -12,14 +12,17 @@
// <cmath>
+#include <array>
#include <cmath>
#include <limits>
#include <type_traits>
#include <cassert>
+#include "fp_compare.h"
#include "test_macros.h"
#include "hexfloat.h"
#include "truncate_fp.h"
+#include "type_algorithms.h"
// convertible to int/float/double/etc
template <class T, int N=0>
@@ -1113,6 +1116,44 @@ void test_fmin()
assert(std::fmin(1,0) == 0);
}
+struct TestHypot3 {
+ template <class Real>
+ static void operator()() {
+ const auto check = [](Real elem, Real abs_tol) {
+ assert(std::isfinite(std::hypot(elem, Real(0), Real(0))));
+ assert(fptest_close(std::hypot(elem, Real(0), Real(0)), elem, abs_tol));
+ assert(std::isfinite(std::hypot(elem, elem, Real(0))));
+ assert(fptest_close(std::hypot(elem, elem, Real(0)), std::sqrt(Real(2)) * elem, abs_tol));
+ assert(std::isfinite(std::hypot(elem, elem, elem)));
+ assert(fptest_close(std::hypot(elem, elem, elem), std::sqrt(Real(3)) * elem, abs_tol));
+ };
+
+ { // check for overflow
+ const auto [elem, abs_tol] = []() -> std::array<Real, 2> {
+ if constexpr (std::is_same_v<Real, float>)
+ return {1e20f, 1e16f};
+ else if constexpr (std::is_same_v<Real, double>)
+ return {1e300, 1e287};
+ else // long double
+ return {1e4000l, 1e3985l};
+ }();
+ check(elem, abs_tol);
+ }
+
+ { // check for underflow
+ const auto [elem, abs_tol] = []() -> std::array<Real, 2> {
+ if constexpr (std::is_same_v<Real, float>)
+ return {1e-20f, 1e-24f};
+ else if constexpr (std::is_same_v<Real, double>)
+ return {1e-287, 1e-300};
+ else // long double
+ return {1e-3985l, 1e-4000l};
+ }();
+ check(elem, abs_tol);
+ }
+ }
+};
+
void test_hypot()
{
static_assert((std::is_same<decltype(std::hypot((float)0, (float)0)), float>::value), "");
@@ -1136,24 +1177,30 @@ void test_hypot()
assert(std::hypot(3,4) == 5);
#if TEST_STD_VER > 14
- static_assert((std::is_same<decltype(std::hypot((float)0, (float)0, (float)0)), float>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (bool)0, (float)0)), double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (unsigned short)0, (double)0)), double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (int)0, (long double)0)), long double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (double)0, (long)0)), double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (long double)0, (unsigned long)0)), long double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (int)0, (long long)0)), double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (int)0, (unsigned long long)0)), double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (double)0, (double)0)), double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (long double)0, (long double)0)), long double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (float)0, (double)0)), double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (float)0, (long double)0)), long double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((float)0, (double)0, (long double)0)), long double>::value), "");
- static_assert((std::is_same<decltype(std::hypot((int)0, (int)0, (int)0)), double>::value), "");
- static_assert((std::is_same<decltype(hypot(Ambiguous(), Ambiguous(), Ambiguous())), Ambiguous>::value), "");
+ // clang-format off
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (float)0, (float)0)), float>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (bool)0, (float)0)), double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (unsigned short)0, (double)0)), double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (int)0, (long double)0)), long double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (double)0, (long)0)), double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (long double)0, (unsigned long)0)), long double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (int)0, (long long)0)), double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (int)0, (unsigned long long)0)), double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (double)0, (double)0)), double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (long double)0, (long double)0)), long double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (float)0, (double)0)), double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (float)0, (long double)0)), long double>));
+ static_assert((std::is_same_v<decltype(std::hypot((float)0, (double)0, (long double)0)), long double>));
+ static_assert((std::is_same_v<decltype(std::hypot((int)0, (int)0, (int)0)), double>));
+ static_assert((std::is_same_v<decltype(hypot(Ambiguous(), Ambiguous(), Ambiguous())), Ambiguous>));
+ // clang-format on
assert(std::hypot(2,3,6) == 7);
assert(std::hypot(1,4,8) == 9);
+
+ // Check for undue over-/underflows of intermediate results.
+ // See discussion at https://github.com/llvm/llvm-project/issues/92782.
+ types::for_each(types::floating_point_types(), TestHypot3());
#endif
}
``````````
</details>
https://github.com/llvm/llvm-project/pull/93350
More information about the libcxx-commits
mailing list