[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