Skip to content

Commit 862c745

Browse files
PaulXiCaoPaulXiCao
committed
implementation+testing
1 parent af06440 commit 862c745

File tree

3 files changed

+124
-24
lines changed

3 files changed

+124
-24
lines changed

‎libcxx/include/__math/hypot.h

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,11 +9,16 @@
99
#ifndef _LIBCPP___MATH_HYPOT_H
1010
#define _LIBCPP___MATH_HYPOT_H
1111

12+
#include <__algorithm/max.h>
1213
#include <__config>
14+
#include <__math/abs.h>
15+
#include <__math/roots.h>
1316
#include <__type_traits/enable_if.h>
1417
#include <__type_traits/is_arithmetic.h>
1518
#include <__type_traits/is_same.h>
1619
#include <__type_traits/promote.h>
20+
#include <array>
21+
#include <limits>
1722

1823
#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
1924
# pragma GCC system_header
@@ -41,6 +46,79 @@ inline _LIBCPP_HIDE_FROM_ABI typename __promote<_A1, _A2>::type hypot(_A1 __x, _
4146
return __math::hypot((__result_type)__x, (__result_type)__y);
4247
}
4348

49+
#if _LIBCPP_STD_VER >= 17
50+
template <class _Real>
51+
struct __hypot_factors {
52+
_Real __threshold;
53+
_Real __scale_xyz;
54+
_Real __scale_M;
55+
};
56+
57+
// returns [underflow_factors, overflow_factors]
58+
template <class _Real>
59+
std::array<__hypot_factors<_Real>, 2> __create_factors() {
60+
static_assert(std::numeric_limits<_Real>::is_iec559);
61+
62+
__hypot_factors<_Real> __underflow, __overflow;
63+
if constexpr (std::is_same_v<_Real, float>) {
64+
static_assert(-125 == std::numeric_limits<_Real>::min_exponent);
65+
static_assert(+128 == std::numeric_limits<_Real>::max_exponent);
66+
__underflow = __hypot_factors<_Real>{0x1.0p-62f, 0x1.0p70f, 0x1.0p-70f};
67+
__overflow = __hypot_factors<_Real>{0x1.0p62f, 0x1.0p-70f, 0x1.0p+70f};
68+
} else if constexpr (std::is_same_v<_Real, double>) {
69+
static_assert(-1021 == std::numeric_limits<_Real>::min_exponent);
70+
static_assert(+1024 == std::numeric_limits<_Real>::max_exponent);
71+
__underflow = __hypot_factors<_Real>{0x1.0p-510, 0x1.0p600, 0x1.0p-600};
72+
__overflow = __hypot_factors<_Real>{0x1.0p510, 0x1.0p-600, 0x1.0p+600};
73+
} else { // long double
74+
static_assert(std::is_same_v<_Real, long double>);
75+
static_assert(-16'381 == std::numeric_limits<_Real>::min_exponent);
76+
static_assert(+16'384 == std::numeric_limits<_Real>::max_exponent);
77+
__underflow = __hypot_factors<_Real>{0x1.0p-8'190l, 0x1.0p9'000l, 0x1.0p-9'000l};
78+
__overflow = __hypot_factors<_Real>{0x1.0p8'190l, 0x1.0p-9'000l, 0x1.0p+9'000l};
79+
}
80+
return {__underflow, __overflow};
81+
}
82+
83+
template <class _Real>
84+
_Real __hypot(_Real __x, _Real __y, _Real __z) {
85+
const auto [__underflow, __overflow] = __create_factors<_Real>();
86+
_Real __M = std::max({__math::fabs(__x), __math::fabs(__y), __math::fabs(__z)});
87+
if (__M > __overflow.__threshold) { // x*x + y*y + z*z might overflow
88+
__x *= __overflow.__scale_xyz;
89+
__y *= __overflow.__scale_xyz;
90+
__z *= __overflow.__scale_xyz;
91+
__M = __overflow.__scale_M;
92+
} else if (__M < __underflow.__threshold) { // x*x + y*y + z*z might underflow
93+
__x *= __underflow.__scale_xyz;
94+
__y *= __underflow.__scale_xyz;
95+
__z *= __underflow.__scale_xyz;
96+
__M = __underflow.__scale_M;
97+
} else
98+
__M = 1;
99+
return __M * __math::sqrt(__x * __x + __y * __y + __z * __z);
100+
}
101+
102+
inline _LIBCPP_HIDE_FROM_ABI float hypot(float __x, float __y, float __z) { return __hypot(__x, __y, __z); }
103+
104+
inline _LIBCPP_HIDE_FROM_ABI double hypot(double __x, double __y, double __z) { return __hypot(__x, __y, __z); }
105+
106+
inline _LIBCPP_HIDE_FROM_ABI long double hypot(long double __x, long double __y, long double __z) {
107+
return __hypot(__x, __y, __z);
108+
}
109+
110+
template <class _A1,
111+
class _A2,
112+
class _A3,
113+
std::enable_if_t< is_arithmetic_v<_A1> && is_arithmetic_v<_A2> && is_arithmetic_v<_A3>, int> = 0 >
114+
_LIBCPP_HIDE_FROM_ABI typename __promote<_A1, _A2, _A3>::type hypot(_A1 __x, _A2 __y, _A3 __z) _NOEXCEPT {
115+
using __result_type = typename __promote<_A1, _A2, _A3>::type;
116+
static_assert(!(
117+
std::is_same_v<_A1, __result_type> && std::is_same_v<_A2, __result_type> && std::is_same_v<_A3, __result_type>));
118+
return __hypot(static_cast<__result_type>(__x), static_cast<__result_type>(__y), static_cast<__result_type>(__z));
119+
}
120+
#endif
121+
44122
} // namespace __math
45123

46124
_LIBCPP_END_NAMESPACE_STD

‎libcxx/include/cmath

Lines changed: 1 addition & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -305,6 +305,7 @@ constexpr long double lerp(long double a, long double b, long double t) noexcept
305305
*/
306306

307307
#include <__config>
308+
#include <__math/hypot.h>
308309
#include <__type_traits/enable_if.h>
309310
#include <__type_traits/is_arithmetic.h>
310311
#include <__type_traits/is_constant_evaluated.h>
@@ -544,30 +545,6 @@ using ::scalbnl _LIBCPP_USING_IF_EXISTS;
544545
using ::tgammal _LIBCPP_USING_IF_EXISTS;
545546
using ::truncl _LIBCPP_USING_IF_EXISTS;
546547

547-
#if _LIBCPP_STD_VER >= 17
548-
inline _LIBCPP_HIDE_FROM_ABI float hypot(float __x, float __y, float __z) {
549-
return sqrt(__x * __x + __y * __y + __z * __z);
550-
}
551-
inline _LIBCPP_HIDE_FROM_ABI double hypot(double __x, double __y, double __z) {
552-
return sqrt(__x * __x + __y * __y + __z * __z);
553-
}
554-
inline _LIBCPP_HIDE_FROM_ABI long double hypot(long double __x, long double __y, long double __z) {
555-
return sqrt(__x * __x + __y * __y + __z * __z);
556-
}
557-
558-
template <class _A1, class _A2, class _A3>
559-
inline _LIBCPP_HIDE_FROM_ABI
560-
typename enable_if_t< is_arithmetic<_A1>::value && is_arithmetic<_A2>::value && is_arithmetic<_A3>::value,
561-
__promote<_A1, _A2, _A3> >::type
562-
hypot(_A1 __lcpp_x, _A2 __lcpp_y, _A3 __lcpp_z) _NOEXCEPT {
563-
typedef typename __promote<_A1, _A2, _A3>::type __result_type;
564-
static_assert(
565-
!(is_same<_A1, __result_type>::value && is_same<_A2, __result_type>::value && is_same<_A3, __result_type>::value),
566-
"");
567-
return std::hypot((__result_type)__lcpp_x, (__result_type)__lcpp_y, (__result_type)__lcpp_z);
568-
}
569-
#endif
570-
571548
template <class _A1, __enable_if_t<is_floating_point<_A1>::value, int> = 0>
572549
_LIBCPP_HIDE_FROM_ABI _LIBCPP_CONSTEXPR bool __constexpr_isnan(_A1 __lcpp_x) _NOEXCEPT {
573550
#if __has_builtin(__builtin_isnan)

‎libcxx/test/std/numerics/c.math/cmath.pass.cpp

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,14 +12,17 @@
1212

1313
// <cmath>
1414

15+
#include <array>
1516
#include <cmath>
1617
#include <limits>
1718
#include <type_traits>
1819
#include <cassert>
1920

21+
#include "fp_compare.h"
2022
#include "test_macros.h"
2123
#include "hexfloat.h"
2224
#include "truncate_fp.h"
25+
#include "type_algorithms.h"
2326

2427
// convertible to int/float/double/etc
2528
template <class T, int N=0>
@@ -1113,6 +1116,44 @@ void test_fmin()
11131116
assert(std::fmin(1,0) == 0);
11141117
}
11151118

1119+
struct TestHypot3 {
1120+
template <class Real>
1121+
static void operator()() {
1122+
const auto check = [](Real elem, Real abs_tol) {
1123+
assert(std::isfinite(std::hypot(elem, Real(0), Real(0))));
1124+
assert(fptest_close(std::hypot(elem, Real(0), Real(0)), elem, abs_tol));
1125+
assert(std::isfinite(std::hypot(elem, elem, Real(0))));
1126+
assert(fptest_close(std::hypot(elem, elem, Real(0)), std::sqrt(Real(2)) * elem, abs_tol));
1127+
assert(std::isfinite(std::hypot(elem, elem, elem)));
1128+
assert(fptest_close(std::hypot(elem, elem, elem), std::sqrt(Real(3)) * elem, abs_tol));
1129+
};
1130+
1131+
{ // check for overflow
1132+
const auto [elem, abs_tol] = []() -> std::array<Real, 2> {
1133+
if constexpr (std::is_same_v<Real, float>)
1134+
return {1e20f, 1e16f};
1135+
else if constexpr (std::is_same_v<Real, double>)
1136+
return {1e300, 1e287};
1137+
else // long double
1138+
return {1e4000l, 1e3985l};
1139+
}();
1140+
check(elem, abs_tol);
1141+
}
1142+
1143+
{ // check for underflow
1144+
const auto [elem, abs_tol] = []() -> std::array<Real, 2> {
1145+
if constexpr (std::is_same_v<Real, float>)
1146+
return {1e-20f, 1e-24f};
1147+
else if constexpr (std::is_same_v<Real, double>)
1148+
return {1e-287, 1e-300};
1149+
else // long double
1150+
return {1e-3985l, 1e-4000l};
1151+
}();
1152+
check(elem, abs_tol);
1153+
}
1154+
}
1155+
};
1156+
11161157
void test_hypot()
11171158
{
11181159
static_assert((std::is_same<decltype(std::hypot((float)0, (float)0)), float>::value), "");
@@ -1156,6 +1197,10 @@ void test_hypot()
11561197

11571198
assert(std::hypot(2,3,6) == 7);
11581199
assert(std::hypot(1,4,8) == 9);
1200+
1201+
// Check for undue over-/underflows of intermediate results.
1202+
// See discussion at https://github.com/llvm/llvm-project/issues/92782.
1203+
types::for_each(types::floating_point_types(), TestHypot3());
11591204
#endif
11601205
}
11611206

0 commit comments

Comments
 (0)