Skip to content

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

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.

Sign up for GitHub

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

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

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3554b85
[libcxx] makes `tanh(complex<float>)` work for large values
cjdb Jan 6, 2025
976a693
adds tests for lowest floating-point values
cjdb Jan 8, 2025
2e9b062
applies clang-format
cjdb Jan 8, 2025
657cfda
finishes clang-format
cjdb Jan 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 30 additions & 6 deletions libcxx/include/complex
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1236,6 +1236,32 @@ _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())) {
Expand All @@ -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.));
Copy link
Contributor

@lntue lntue Jan 9, 2025

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.

return complex<_Tp>(__2rsh / __d, std::__sin2(__x.imag()) / __d);
}

// asin
Expand Down
Loading
Loading