Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[libcxx] makes tanh(complex<float>) work for large values #122194

Open
wants to merge 4 commits into
base: main
Choose a base branch
from

Conversation

cjdb
Copy link
Contributor

@cjdb cjdb commented Jan 8, 2025

By attempting to compute sin(2x) and cos(2x), we were inadvertently producing complex<float>(NaN, NaN) for any |x| > 2^127, since 2x would be computed as inf, and sin(inf) == cos(inf) == NaN.

This commit relies on trig identities to sidestep this issue and always produce a valid answer for large values of tanh. The test cases only handled double, so they have been expanded to test float and long double as well.

cjdb added 2 commits January 8, 2025 23:45
By attempting to compute `sin(2x)` and `cos(2x)`, we were inadvertently
producing `complex<float>(NaN, NaN)` for any `|x| > 2^127`, since `2x`
would be computed as `inf`, and `sin(inf) == cos(inf) == NaN`.

This commit relies on trig identities to sidestep this issue and always
produce a valid answer for large values of `tanh`. The test cases only
handled `double`, so they have been expanded to test `float` and `long
double` as well.
@cjdb cjdb added libc++ libc++ C++ Standard Library. Not GNU libstdc++. Not libc++abi. floating-point Floating-point math labels Jan 8, 2025
@cjdb cjdb requested a review from lntue January 8, 2025 23:48
@cjdb cjdb requested a review from a team as a code owner January 8, 2025 23:48
@llvmbot
Copy link
Member

llvmbot commented Jan 8, 2025

@llvm/pr-subscribers-libcxx

Author: Christopher Di Bella (cjdb)

Changes

By attempting to compute sin(2x) and cos(2x), we were inadvertently producing complex&lt;float&gt;(NaN, NaN) for any |x| &gt; 2^127, since 2x would be computed as inf, and sin(inf) == cos(inf) == NaN.

This commit relies on trig identities to sidestep this issue and always produce a valid answer for large values of tanh. The test cases only handled double, so they have been expanded to test float and long double as well.


Patch is 117.15 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/122194.diff

28 Files Affected:

  • (modified) libcxx/include/complex (+34-6)
  • (modified) libcxx/test/std/numerics/complex.number/cases.h (+198-172)
  • (modified) libcxx/test/std/numerics/complex.number/complex.ops/complex_divide_complex.pass.cpp (+31-6)
  • (modified) libcxx/test/std/numerics/complex.number/complex.ops/complex_times_complex.pass.cpp (+26-6)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/acos.pass.cpp (+39-29)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp (+33-30)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/asin.pass.cpp (+29-26)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp (+32-29)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/atan.pass.cpp (+9-6)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/atanh.pass.cpp (+34-31)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/cos.pass.cpp (+8-5)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/cosh.pass.cpp (+29-19)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/exp.pass.cpp (+24-21)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/log.pass.cpp (+27-24)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/log10.pass.cpp (+7-4)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/pow_complex_complex.pass.cpp (+7-4)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/pow_complex_scalar.pass.cpp (+7-4)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/pow_scalar_complex.pass.cpp (+7-4)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/sin.pass.cpp (+9-6)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/sinh.pass.cpp (+32-22)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/sqrt.pass.cpp (+53-19)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/tan.pass.cpp (+9-6)
  • (modified) libcxx/test/std/numerics/complex.number/complex.transcendentals/tanh.pass.cpp (+28-19)
  • (modified) libcxx/test/std/numerics/complex.number/complex.value.ops/abs.pass.cpp (+21-4)
  • (modified) libcxx/test/std/numerics/complex.number/complex.value.ops/arg.pass.cpp (+34-29)
  • (modified) libcxx/test/std/numerics/complex.number/complex.value.ops/norm.pass.cpp (+9-4)
  • (modified) libcxx/test/std/numerics/complex.number/complex.value.ops/polar.pass.cpp (+11-6)
  • (modified) libcxx/test/std/numerics/complex.number/complex.value.ops/proj.pass.cpp (+14-9)
diff --git a/libcxx/include/complex b/libcxx/include/complex
index df18159595b34d..5f425206100ce2 100644
--- a/libcxx/include/complex
+++ b/libcxx/include/complex
@@ -1236,6 +1236,36 @@ _LIBCPP_HIDE_FROM_ABI complex<_Tp> cosh(const complex<_Tp>& __x) {
 
 // tanh
 
+template<class _Tp>
+_LIBCPP_HIDE_FROM_ABI _Tp __sin2(const _Tp __x) noexcept
+{
+  static_assert(std::is_arithmetic<_Tp>::value, "requires an arithmetic type");
+  return 2 * std::sin(__x) * std::cos(__x);
+}
+
+template<class _Tp>
+_LIBCPP_HIDE_FROM_ABI _Tp __sinh2(const _Tp __x) noexcept
+{
+  static_assert(std::is_arithmetic<_Tp>::value, "requires an arithmetic type");
+  return 2 * std::sinh(__x) * std::cosh(__x);
+}
+
+template<class _Tp>
+_LIBCPP_HIDE_FROM_ABI _Tp __cos2(const _Tp __x) noexcept
+{
+  static_assert(std::is_arithmetic<_Tp>::value, "requires an arithmetic type");
+  const _Tp __cos = std::cos(__x);
+  return 2 * __cos * __cos - 1;
+}
+
+template<class _Tp>
+_LIBCPP_HIDE_FROM_ABI _Tp __cosh2(const _Tp __x) noexcept
+{
+  static_assert(std::is_arithmetic<_Tp>::value, "requires an arithmetic type");
+  const _Tp __cosh = std::cosh(__x);
+  return 2 * __cosh * __cosh - 1;
+}
+
 template <class _Tp>
 _LIBCPP_HIDE_FROM_ABI complex<_Tp> tanh(const complex<_Tp>& __x) {
   if (std::isinf(__x.real())) {
@@ -1245,13 +1275,11 @@ _LIBCPP_HIDE_FROM_ABI complex<_Tp> tanh(const complex<_Tp>& __x) {
   }
   if (std::isnan(__x.real()) && __x.imag() == 0)
     return __x;
-  _Tp __2r(_Tp(2) * __x.real());
-  _Tp __2i(_Tp(2) * __x.imag());
-  _Tp __d(std::cosh(__2r) + std::cos(__2i));
-  _Tp __2rsh(std::sinh(__2r));
+  _Tp __d(std::__cosh2(__x.real()) + std::__cos2(__x.imag()));
+  _Tp __2rsh(std::__sinh2(__x.real()));
   if (std::isinf(__2rsh) && std::isinf(__d))
-    return complex<_Tp>(__2rsh > _Tp(0) ? _Tp(1) : _Tp(-1), __2i > _Tp(0) ? _Tp(0) : _Tp(-0.));
-  return complex<_Tp>(__2rsh / __d, std::sin(__2i) / __d);
+    return complex<_Tp>(__2rsh > _Tp(0) ? _Tp(1) : _Tp(-1), __x.imag() > _Tp(0) ? _Tp(0) : _Tp(-0.));
+  return complex<_Tp>(__2rsh / __d, std::__sin2(__x.imag()) / __d);
 }
 
 // asin
diff --git a/libcxx/test/std/numerics/complex.number/cases.h b/libcxx/test/std/numerics/complex.number/cases.h
index b360e1423ff576..37d327450311ae 100644
--- a/libcxx/test/std/numerics/complex.number/cases.h
+++ b/libcxx/test/std/numerics/complex.number/cases.h
@@ -15,184 +15,201 @@
 
 #include <cassert>
 #include <complex>
+#include <limits>
 #include <type_traits>
 
 #include "test_macros.h"
 
-TEST_CONSTEXPR_CXX20 const std::complex<double> testcases[] =
+template<class T>
+TEST_CONSTEXPR_CXX20 const std::complex<T> testcases[] =
 {
-    std::complex<double>( 1.e-6,  1.e-6),
-    std::complex<double>(-1.e-6,  1.e-6),
-    std::complex<double>(-1.e-6, -1.e-6),
-    std::complex<double>( 1.e-6, -1.e-6),
-
-    std::complex<double>( 1.e+6,  1.e-6),
-    std::complex<double>(-1.e+6,  1.e-6),
-    std::complex<double>(-1.e+6, -1.e-6),
-    std::complex<double>( 1.e+6, -1.e-6),
-
-    std::complex<double>( 1.e-6,  1.e+6),
-    std::complex<double>(-1.e-6,  1.e+6),
-    std::complex<double>(-1.e-6, -1.e+6),
-    std::complex<double>( 1.e-6, -1.e+6),
-
-    std::complex<double>( 1.e+6,  1.e+6),
-    std::complex<double>(-1.e+6,  1.e+6),
-    std::complex<double>(-1.e+6, -1.e+6),
-    std::complex<double>( 1.e+6, -1.e+6),
-
-    std::complex<double>(-0, -1.e-6),
-    std::complex<double>(-0,  1.e-6),
-    std::complex<double>(-0,  1.e+6),
-    std::complex<double>(-0, -1.e+6),
-    std::complex<double>( 0, -1.e-6),
-    std::complex<double>( 0,  1.e-6),
-    std::complex<double>( 0,  1.e+6),
-    std::complex<double>( 0, -1.e+6),
-
-    std::complex<double>(-1.e-6, -0),
-    std::complex<double>( 1.e-6, -0),
-    std::complex<double>( 1.e+6, -0),
-    std::complex<double>(-1.e+6, -0),
-    std::complex<double>(-1.e-6,  0),
-    std::complex<double>( 1.e-6,  0),
-    std::complex<double>( 1.e+6,  0),
-    std::complex<double>(-1.e+6,  0),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(-2, std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(-1, std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(-0.5, std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(-0., std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(+0., std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(0.5, std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(1, std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(2, std::numeric_limits<double>::quiet_NaN()),
-    std::complex<double>(std::numeric_limits<double>::infinity(), std::numeric_limits<double>::quiet_NaN()),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), -std::numeric_limits<double>::infinity()),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity()),
-    std::complex<double>(-2, -std::numeric_limits<double>::infinity()),
-    std::complex<double>(-1, -std::numeric_limits<double>::infinity()),
-    std::complex<double>(-0.5, -std::numeric_limits<double>::infinity()),
-    std::complex<double>(-0., -std::numeric_limits<double>::infinity()),
-    std::complex<double>(+0., -std::numeric_limits<double>::infinity()),
-    std::complex<double>(0.5, -std::numeric_limits<double>::infinity()),
-    std::complex<double>(1, -std::numeric_limits<double>::infinity()),
-    std::complex<double>(2, -std::numeric_limits<double>::infinity()),
-    std::complex<double>(std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity()),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), -2),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), -2),
-    std::complex<double>(-2, -2),
-    std::complex<double>(-1, -2),
-    std::complex<double>(-0.5, -2),
-    std::complex<double>(-0., -2),
-    std::complex<double>(+0., -2),
-    std::complex<double>(0.5, -2),
-    std::complex<double>(1, -2),
-    std::complex<double>(2, -2),
-    std::complex<double>(std::numeric_limits<double>::infinity(), -2),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), -1),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), -1),
-    std::complex<double>(-2, -1),
-    std::complex<double>(-1, -1),
-    std::complex<double>(-0.5, -1),
-    std::complex<double>(-0., -1),
-    std::complex<double>(+0., -1),
-    std::complex<double>(0.5, -1),
-    std::complex<double>(1, -1),
-    std::complex<double>(2, -1),
-    std::complex<double>(std::numeric_limits<double>::infinity(), -1),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), -0.5),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), -0.5),
-    std::complex<double>(-2, -0.5),
-    std::complex<double>(-1, -0.5),
-    std::complex<double>(-0.5, -0.5),
-    std::complex<double>(-0., -0.5),
-    std::complex<double>(+0., -0.5),
-    std::complex<double>(0.5, -0.5),
-    std::complex<double>(1, -0.5),
-    std::complex<double>(2, -0.5),
-    std::complex<double>(std::numeric_limits<double>::infinity(), -0.5),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), -0.),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), -0.),
-    std::complex<double>(-2, -0.),
-    std::complex<double>(-1, -0.),
-    std::complex<double>(-0.5, -0.),
-    std::complex<double>(-0., -0.),
-    std::complex<double>(+0., -0.),
-    std::complex<double>(0.5, -0.),
-    std::complex<double>(1, -0.),
-    std::complex<double>(2, -0.),
-    std::complex<double>(std::numeric_limits<double>::infinity(), -0.),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), +0.),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), +0.),
-    std::complex<double>(-2, +0.),
-    std::complex<double>(-1, +0.),
-    std::complex<double>(-0.5, +0.),
-    std::complex<double>(-0., +0.),
-    std::complex<double>(+0., +0.),
-    std::complex<double>(0.5, +0.),
-    std::complex<double>(1, +0.),
-    std::complex<double>(2, +0.),
-    std::complex<double>(std::numeric_limits<double>::infinity(), +0.),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), 0.5),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), 0.5),
-    std::complex<double>(-2, 0.5),
-    std::complex<double>(-1, 0.5),
-    std::complex<double>(-0.5, 0.5),
-    std::complex<double>(-0., 0.5),
-    std::complex<double>(+0., 0.5),
-    std::complex<double>(0.5, 0.5),
-    std::complex<double>(1, 0.5),
-    std::complex<double>(2, 0.5),
-    std::complex<double>(std::numeric_limits<double>::infinity(), 0.5),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), 1),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), 1),
-    std::complex<double>(-2, 1),
-    std::complex<double>(-1, 1),
-    std::complex<double>(-0.5, 1),
-    std::complex<double>(-0., 1),
-    std::complex<double>(+0., 1),
-    std::complex<double>(0.5, 1),
-    std::complex<double>(1, 1),
-    std::complex<double>(2, 1),
-    std::complex<double>(std::numeric_limits<double>::infinity(), 1),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), 2),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), 2),
-    std::complex<double>(-2, 2),
-    std::complex<double>(-1, 2),
-    std::complex<double>(-0.5, 2),
-    std::complex<double>(-0., 2),
-    std::complex<double>(+0., 2),
-    std::complex<double>(0.5, 2),
-    std::complex<double>(1, 2),
-    std::complex<double>(2, 2),
-    std::complex<double>(std::numeric_limits<double>::infinity(), 2),
-
-    std::complex<double>(std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::infinity()),
-    std::complex<double>(-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity()),
-    std::complex<double>(-2, std::numeric_limits<double>::infinity()),
-    std::complex<double>(-1, std::numeric_limits<double>::infinity()),
-    std::complex<double>(-0.5, std::numeric_limits<double>::infinity()),
-    std::complex<double>(-0., std::numeric_limits<double>::infinity()),
-    std::complex<double>(+0., std::numeric_limits<double>::infinity()),
-    std::complex<double>(0.5, std::numeric_limits<double>::infinity()),
-    std::complex<double>(1, std::numeric_limits<double>::infinity()),
-    std::complex<double>(2, std::numeric_limits<double>::infinity()),
-    std::complex<double>(std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity())
+    std::complex<T>( 1.e-6,  1.e-6),
+    std::complex<T>(-1.e-6,  1.e-6),
+    std::complex<T>(-1.e-6, -1.e-6),
+    std::complex<T>( 1.e-6, -1.e-6),
+
+    std::complex<T>( 1.e+6,  1.e-6),
+    std::complex<T>(-1.e+6,  1.e-6),
+    std::complex<T>(-1.e+6, -1.e-6),
+    std::complex<T>( 1.e+6, -1.e-6),
+
+    std::complex<T>( 1.e-6,  1.e+6),
+    std::complex<T>(-1.e-6,  1.e+6),
+    std::complex<T>(-1.e-6, -1.e+6),
+    std::complex<T>( 1.e-6, -1.e+6),
+
+    std::complex<T>( 1.e+6,  1.e+6),
+    std::complex<T>(-1.e+6,  1.e+6),
+    std::complex<T>(-1.e+6, -1.e+6),
+    std::complex<T>( 1.e+6, -1.e+6),
+
+    std::complex<T>(-0, -1.e-6),
+    std::complex<T>(-0,  1.e-6),
+    std::complex<T>(-0,  1.e+6),
+    std::complex<T>(-0, -1.e+6),
+    std::complex<T>( 0, -1.e-6),
+    std::complex<T>( 0,  1.e-6),
+    std::complex<T>( 0,  1.e+6),
+    std::complex<T>( 0, -1.e+6),
+
+    std::complex<T>(-1.e-6, -0),
+    std::complex<T>( 1.e-6, -0),
+    std::complex<T>( 1.e+6, -0),
+    std::complex<T>(-1.e+6, -0),
+    std::complex<T>(-1.e-6,  0),
+    std::complex<T>( 1.e-6,  0),
+    std::complex<T>( 1.e+6,  0),
+    std::complex<T>(-1.e+6,  0),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(-2, std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(-1, std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(-0.5, std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(-0., std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(+0., std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(0.5, std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(1, std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(2, std::numeric_limits<T>::quiet_NaN()),
+    std::complex<T>(std::numeric_limits<T>::infinity(), std::numeric_limits<T>::quiet_NaN()),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), -std::numeric_limits<T>::infinity()),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity()),
+    std::complex<T>(-2, -std::numeric_limits<T>::infinity()),
+    std::complex<T>(-1, -std::numeric_limits<T>::infinity()),
+    std::complex<T>(-0.5, -std::numeric_limits<T>::infinity()),
+    std::complex<T>(-0., -std::numeric_limits<T>::infinity()),
+    std::complex<T>(+0., -std::numeric_limits<T>::infinity()),
+    std::complex<T>(0.5, -std::numeric_limits<T>::infinity()),
+    std::complex<T>(1, -std::numeric_limits<T>::infinity()),
+    std::complex<T>(2, -std::numeric_limits<T>::infinity()),
+    std::complex<T>(std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity()),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), -2),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), -2),
+    std::complex<T>(-2, -2),
+    std::complex<T>(-1, -2),
+    std::complex<T>(-0.5, -2),
+    std::complex<T>(-0., -2),
+    std::complex<T>(+0., -2),
+    std::complex<T>(0.5, -2),
+    std::complex<T>(1, -2),
+    std::complex<T>(2, -2),
+    std::complex<T>(std::numeric_limits<T>::infinity(), -2),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), -1),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), -1),
+    std::complex<T>(-2, -1),
+    std::complex<T>(-1, -1),
+    std::complex<T>(-0.5, -1),
+    std::complex<T>(-0., -1),
+    std::complex<T>(+0., -1),
+    std::complex<T>(0.5, -1),
+    std::complex<T>(1, -1),
+    std::complex<T>(2, -1),
+    std::complex<T>(std::numeric_limits<T>::infinity(), -1),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), -0.5),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), -0.5),
+    std::complex<T>(-2, -0.5),
+    std::complex<T>(-1, -0.5),
+    std::complex<T>(-0.5, -0.5),
+    std::complex<T>(-0., -0.5),
+    std::complex<T>(+0., -0.5),
+    std::complex<T>(0.5, -0.5),
+    std::complex<T>(1, -0.5),
+    std::complex<T>(2, -0.5),
+    std::complex<T>(std::numeric_limits<T>::infinity(), -0.5),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), -0.),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), -0.),
+    std::complex<T>(-2, -0.),
+    std::complex<T>(-1, -0.),
+    std::complex<T>(-0.5, -0.),
+    std::complex<T>(-0., -0.),
+    std::complex<T>(+0., -0.),
+    std::complex<T>(0.5, -0.),
+    std::complex<T>(1, -0.),
+    std::complex<T>(2, -0.),
+    std::complex<T>(std::numeric_limits<T>::infinity(), -0.),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), +0.),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), +0.),
+    std::complex<T>(-2, +0.),
+    std::complex<T>(-1, +0.),
+    std::complex<T>(-0.5, +0.),
+    std::complex<T>(-0., +0.),
+    std::complex<T>(+0., +0.),
+    std::complex<T>(0.5, +0.),
+    std::complex<T>(1, +0.),
+    std::complex<T>(2, +0.),
+    std::complex<T>(std::numeric_limits<T>::infinity(), +0.),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), 0.5),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), 0.5),
+    std::complex<T>(-2, 0.5),
+    std::complex<T>(-1, 0.5),
+    std::complex<T>(-0.5, 0.5),
+    std::complex<T>(-0., 0.5),
+    std::complex<T>(+0., 0.5),
+    std::complex<T>(0.5, 0.5),
+    std::complex<T>(1, 0.5),
+    std::complex<T>(2, 0.5),
+    std::complex<T>(std::numeric_limits<T>::infinity(), 0.5),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), 1),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), 1),
+    std::complex<T>(-2, 1),
+    std::complex<T>(-1, 1),
+    std::complex<T>(-0.5, 1),
+    std::complex<T>(-0., 1),
+    std::complex<T>(+0., 1),
+    std::complex<T>(0.5, 1),
+    std::complex<T>(1, 1),
+    std::complex<T>(2, 1),
+    std::complex<T>(std::numeric_limits<T>::infinity(), 1),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), 2),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), 2),
+    std::complex<T>(-2, 2),
+    std::complex<T>(-1, 2),
+    std::complex<T>(-0.5, 2),
+    std::complex<T>(-0., 2),
+    std::complex<T>(+0., 2),
+    std::complex<T>(0.5, 2),
+    std::complex<T>(1, 2),
+    std::complex<T>(2, 2),
+    std::complex<T>(std::numeric_limits<T>::infinity(), 2),
+
+    std::complex<T>(std::numeric_limits<T>::quiet_NaN(), std::numeric_limits<T>::infinity()),
+    std::complex<T>(-std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity()),
+    std::complex<T>(-2, std::numeric_limits<T>::infinity()),
+    std::complex<T>(-1, std::numeric_limits<T>::infinity()),
+    std::complex<T>(-0.5, std::numeric_limits<T>::infinity()),
+    std::complex<T>(-0., std::numeric_limits<T>::infinity()),
+    std::complex<T>(+0., std::numeric_limits<T>::infinity()),
+    std::complex<T>(0.5, std::numeric_limits<T>::infinity()),
+    std::complex<T>(1, std::numeric_limits<T>::infinity()),
+    std::complex<T>(2, std::numeric_limits<T>::infinity()),
+    std::complex<T>(std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity()),
+
+    std::complex<T>(std::numeric_limits<T>::max(), 1),
+    std::complex<T>(std::numeric_limits<T>::max(), -1),
+    std::complex<T>(std::numeric_limits<T>::lowest(), 1),
+    std::complex<T>(std::numeric_limits<T>::lowest(), -1),
+
+    std::complex<T>(1, std::numeric_limits<T>::max()),
+    std::complex<T>(1, std::numeric_limits<T>::lowest()),
+    std::complex<T>(-1, std::numeric_limits<T>::max()),
+    std::complex<T>(-1, std::numeric_limits<T>::lowest()),
+
+    std::complex<T>(std::numeric_limits<T>::max(), std::numeric_limits<T>::max()),
+    std::complex<T>(std::numeric_limits<T>::max(), std::numeric_limits<T>::lowest()),
+    std::complex<T>(std::numeric_limits<T>::lowest(), std::numeric_limits<T>::max()),
+    std::complex<T>(std::numeric_limits<T>::lowest(), std::numeric_limits<T>::lowest()),
 };
 
-enum {zero, non_zero, inf, NaN, non_zero_nan};
+enum {zero, non_zero, lowest_value, maximum_value, inf, NaN, non_zero_nan};
 
 template <class T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
 TEST_CONSTEXPR_CXX20 bool test_isinf(T v) {
@@ -227,12 +244,17 @@ classify(const std::complex<T>& x)
             return NaN;
         return non_zero_nan;
     }
+    if (x.real() == std::numeric_limits<T>::max() || x.imag() == std::numeric_limits<T>::max())
+        return maximum_value;
+    if (x.real() == std::numeric_limits<T>::lowest() || x.imag() == std::numeric_limits<T>::lowest())
+        return lowest_value;
     return non_zero;
 }
 
+template<class T>
 inline
 int
-classify(double x)
+classify(T x)
 {
     if (x == 0)
         return zero;
@@ -240,6 +262,10 @@ classify(double x)
         return inf;
     if (std::isnan(x))
         return NaN;
+    if (x == std::numeric_limits<T>::max())
+        return maximum_value;
+    if (x == std::numeric_limits<T>::lowest())
+        return lowest_value;
     return non_zero;
 }
 
diff --git a/libcxx/test/std/numerics/complex.number/complex.ops/complex_divide_complex.pass.cpp b/libcxx/test/std/numerics/complex.number/complex.ops/complex_divide_complex.pass.cpp
index d12dfd994b0ae9..bc0bdd0ec261dd 100644
--- a/libcxx/test/std/numerics/complex.number/complex.ops/complex_divide_complex.pass.cpp
+++ b/libcxx/test/std/numerics/complex.number/complex.ops/complex_divide_complex.pass.cpp
@@ -12,7 +12,7 @@
 //   complex<T>
 //   operator/(const complex<T>& lhs, const complex<T>& rhs); // constexpr in C++20
 
-// ADDITIONAL_COMPILE_FLAGS(has-fconstexpr-steps): -fconstexpr-steps=2000000
+// ADDITIONAL_COMPILE_FLAGS(has-fconstexpr-steps): -fconstexpr-steps=2131685
 
 #include <cassert>
 #include <complex>
@@ -34,24 +34,30 @@ test()
     return true;
 }
 
+...
[truncated]

Copy link

github-actions bot commented Jan 8, 2025

✅ With the latest revision this PR passed the C/C++ code formatter.

if (std::isinf(__2rsh) && std::isinf(__d))
return complex<_Tp>(__2rsh > _Tp(0) ? _Tp(1) : _Tp(-1), __2i > _Tp(0) ? _Tp(0) : _Tp(-0.));
return complex<_Tp>(__2rsh / __d, std::sin(__2i) / __d);
return complex<_Tp>(__2rsh > _Tp(0) ? _Tp(1) : _Tp(-1), __x.imag() > _Tp(0) ? _Tp(0) : _Tp(-0.));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of switching to a different formula, which is actually less accurate for the normal range in the current form, I think it is better to do the cutoff a lot earlier. This return is actually correct (for the default rounding mode) when:

$$e^{\left|\_\_2r\right|} > \frac{4}{\text{std::numeric\_limits<\_Tp>::epsilon()}}.$$

so for the real part, to prevent overflow, some simple check like:

  std::fabs(__x.real()) >= std::numeric_limits<_Tp>::digits

would be sufficient.

Now for the imaginary part, when

   std::fabs(__x.imag()) > std::numeric_limits<_Tp>::max() / 2

we will need to expand double-angle formulas to compute in __x.imag() only. In that case, I would use the following formula to reduce the cancellation:

$$\begin{align*} \tanh(re + i \cdot im) &= \frac{\sinh(2 \cdot re)}{\cosh(2 \cdot re) + \cos(2 \cdot im) } + i \cdot \frac{\sin(2 \cdot im)}{\cosh(2 \cdot re) + \cos(2 \cdot im) } \\\ &= \frac{0.5 \cdot \sinh(2 \cdot re)}{\sinh^2(re) + \cos^2(im)} + i \cdot \frac{\sin(im) \cdot \cos(im)}{\sinh^2(re) + \cos^2(im) } \end{align*}$$

And use the original formula/implementation for the default case.

Copy link
Contributor

@philnik777 philnik777 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think an approach similar to #99677 would be better. That way we avoid a bunch of code and most likely get a better implementation.

@cjdb
Copy link
Contributor Author

cjdb commented Jan 9, 2025

I think an approach similar to #99677 would be better. That way we avoid a bunch of code and most likely get a better implementation.

When will #99677 land? I'd be happy to reduce this PR to just changes to the tests, if it does.

@philnik777
Copy link
Contributor

I think an approach similar to #99677 would be better. That way we avoid a bunch of code and most likely get a better implementation.

When will #99677 land? I'd be happy to reduce this PR to just changes to the tests, if it does.

I don't know when it'll land, but I'm also not touching tanh in there right now. I'm only talking about the approach (that being just calling the libc versions instead of implementing them on our own), not the patch itself.

@cjdb
Copy link
Contributor Author

cjdb commented Jan 9, 2025

I think an approach similar to #99677 would be better. That way we avoid a bunch of code and most likely get a better implementation.

When will #99677 land? I'd be happy to reduce this PR to just changes to the tests, if it does.

I don't know when it'll land, but I'm also not touching tanh in there right now. I'm only talking about the approach (that being just calling the libc versions instead of implementing them on our own), not the patch itself.

Ah, you're editing atanh, which I misread as tanh. In that case, yes, I'll try that out (@lntue I think this answers your question about using the C99 complex functions).

@lntue
Copy link
Contributor

lntue commented Jan 9, 2025

I think an approach similar to #99677 would be better. That way we avoid a bunch of code and most likely get a better implementation.

When will #99677 land? I'd be happy to reduce this PR to just changes to the tests, if it does.

I don't know when it'll land, but I'm also not touching tanh in there right now. I'm only talking about the approach (that being just calling the libc versions instead of implementing them on our own), not the patch itself.

Ah, you're editing atanh, which I misread as tanh. In that case, yes, I'll try that out (@lntue I think this answers your question about using the C99 complex functions).

So I assume we can do the same for std::sqrt(std::complex<T>) --> csqrt*()?

@cjdb
Copy link
Contributor Author

cjdb commented Jan 9, 2025

Hmm, it seems __builtin_tanhf produces different results to what we do today, and maintaining backwards compatibility complicates the implementation quite a bit (relative to #99677).

I can delete that complexity and adjust the tests, but it won't be transparent, thanks to the clang-format changes.

@cjdb
Copy link
Contributor Author

cjdb commented Jan 9, 2025

I think an approach similar to #99677 would be better. That way we avoid a bunch of code and most likely get a better implementation.

When will #99677 land? I'd be happy to reduce this PR to just changes to the tests, if it does.

I don't know when it'll land, but I'm also not touching tanh in there right now. I'm only talking about the approach (that being just calling the libc versions instead of implementing them on our own), not the patch itself.

Ah, you're editing atanh, which I misread as tanh. In that case, yes, I'll try that out (@lntue I think this answers your question about using the C99 complex functions).

So I assume we can do the same for std::sqrt(std::complex<T>) --> csqrt*()?

Yes, I'll do that in a separate PR.

@cjdb
Copy link
Contributor Author

cjdb commented Jan 9, 2025

Currently blocked by #122391.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
floating-point Floating-point math libc++ libc++ C++ Standard Library. Not GNU libstdc++. Not libc++abi.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants