[libcxx-commits] [libcxx] [libc++][math] Fix undue overflowing of `std::hypot(x, y, z)` (PR #93350)

[Louis Dionne via libcxx-commits]

================
@@ -41,6 +46,83 @@ 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>);
+    if constexpr (sizeof(_Real) == sizeof(double))
+      return static_cast<std::array<__hypot_factors<_Real>, 2>>(__create_factors<double>());
+    else {
+      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); }
----------------
ldionne wrote:

Those qualifications are to prevent ADL from kicking in, and also for consistency (since ADL doesn't really do anything in this case).

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