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

via libcxx-commits libcxx-commits at lists.llvm.org
Thu Jan 9 08:22:21 PST 2025


================
@@ -1245,13 +1271,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.));
----------------
lntue wrote:

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:
```math
  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:
```math
\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.

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


More information about the libcxx-commits mailing list