[libcxx] Improve accuracy of complex asinh and acosh

[miyuki]

Summary

Currently std::asinh and std::acosh use std::pow to compute x^2. This
results in a significant error when computing e.g. asinh(i) or
acosh(-1).

This patch expresses x^2 directly via x.real() and x.imag(), like it
is done in libstdc++/glibc, and adds tests that checks the accuracy.

[miyuki]

test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp

59

Ideally, I would prefer some approximate comparison instead of == here, but is_about is not suitable for arguments that have zero or near-zero sum. That's why I used exact comparison (several other cases in this file do the same). This test passes on x86_64, so hopefully it will work the same way in any IEEE-compliant environment.

[miyuki]

ping

[miyuki]

ping^2

[miyuki]

ping^3

[mclow.lists]

This all looks good to me.
I think that one more test should be added - and that's one that tests __sqr directly.
Since that's not a public routine, the test should go in "test/libcxx/numerics/complex.number"

[miyuki]

Added a test for __sqr

[miyuki]

ping

[mclow.lists]

LGTM. Do you need me to commit this?

[miyuki]

No, I'll commit it myself. Thank you for the review.

include/complex

<
     complex<typename __promote<_Tp, _Up>::type>
 >::type
 pow(const _Tp& __x, const complex<_Up>& __y)
 {
     typedef complex<typename __promote<_Tp, _Up>::type> result_type;
     return _VSTD::pow(result_type(__x), result_type(__y));
 }
 
+// __sqr, computes pow(x, 2)
+
+template<class _Tp>
+inline _LIBCPP_INLINE_VISIBILITY
+complex<_Tp>
+__sqr(const complex<_Tp>& __x)
+{
+    return complex<_Tp>((__x.real() - __x.imag()) * (__x.real() + __x.imag()),
+                        _Tp(2) * __x.real() * __x.imag());
+}
+
 // asinh
 
 template<class _Tp>
 complex<_Tp>
 asinh(const complex<_Tp>& __x)
 {
     const _Tp __pi(atan2(+0., -0.));
     if (__libcpp_isinf_or_builtin(__x.real()))
if (__libcpp_isnan_or_builtin(__x.real()))
         if (__libcpp_isinf_or_builtin(__x.imag()))
             return complex<_Tp>(__x.imag(), __x.real());
         if (__x.imag() == 0)
             return __x;
         return complex<_Tp>(__x.real(), __x.real());
     }
     if (__libcpp_isinf_or_builtin(__x.imag()))
         return complex<_Tp>(copysign(__x.imag(), __x.real()), copysign(__pi/_Tp(2), __x.imag()));
-    complex<_Tp> __z = log(__x + sqrt(pow(__x, _Tp(2)) + _Tp(1)));
+    complex<_Tp> __z = log(__x + sqrt(__sqr(__x) + _Tp(1)));
     return complex<_Tp>(copysign(__z.real(), __x.real()), copysign(__z.imag(), __x.imag()));
 }
 
 // acosh
 
 template<class _Tp>
 complex<_Tp>
 acosh(const complex<_Tp>& __x)
acosh(const complex<_Tp>& __x)
     if (__libcpp_isnan_or_builtin(__x.real()))
     {
         if (__libcpp_isinf_or_builtin(__x.imag()))
             return complex<_Tp>(abs(__x.imag()), __x.real());
         return complex<_Tp>(__x.real(), __x.real());
     }
     if (__libcpp_isinf_or_builtin(__x.imag()))
         return complex<_Tp>(abs(__x.imag()), copysign(__pi/_Tp(2), __x.imag()));
-    complex<_Tp> __z = log(__x + sqrt(pow(__x, _Tp(2)) - _Tp(1)));
+    complex<_Tp> __z = log(__x + sqrt(__sqr(__x) - _Tp(1)));
     return complex<_Tp>(copysign(__z.real(), _Tp(0)), copysign(__z.imag(), __x.imag()));
 }
 
 // atanh
 
 template<class _Tp>
 complex<_Tp>
 atanh(const complex<_Tp>& __x)
if (__libcpp_isnan_or_builtin(__x.real()))
         if (__libcpp_isinf_or_builtin(__x.imag()))
             return complex<_Tp>(__x.real(), -__x.imag());
         return complex<_Tp>(__x.real(), __x.real());
     }
     if (__libcpp_isinf_or_builtin(__x.imag()))
         return complex<_Tp>(__pi/_Tp(2), -__x.imag());
     if (__x.real() == 0 && (__x.imag() == 0 || isnan(__x.imag())))
         return complex<_Tp>(__pi/_Tp(2), -__x.imag());
-    complex<_Tp> __z = log(__x + sqrt(pow(__x, _Tp(2)) - _Tp(1)));
+    complex<_Tp> __z = log(__x + sqrt(__sqr(__x) - _Tp(1)));
     if (signbit(__x.imag()))
         return complex<_Tp>(abs(__z.imag()), abs(__z.real()));
     return complex<_Tp>(abs(__z.imag()), -abs(__z.real()));
 }
 
 // atan
 
 template<class _Tp>

test/libcxx/numerics/complex.number/__sqr.pass.cpp

+//===----------------------------------------------------------------------===//
+//
+//                     The LLVM Compiler Infrastructure
+//
+// This file is dual licensed under the MIT and the University of Illinois Open
+// Source Licenses. See LICENSE.TXT for details.
+//
+//===----------------------------------------------------------------------===//
+
+// <complex>
+
+// template<class T>
+//   complex<T>
+//   __sqr(const complex<T>& x);
+
+#include <complex>
+#include <cassert>
+
+template <class T>
+void
+test()
+{
+    const T tolerance = std::is_same<T, float>::value ? 1.e-6 : 1.e-14;
+
+    typedef std::complex<T> cplx;
+    struct test_case
+    {
+        cplx value;
+        cplx expected;
+    };
+
+    const test_case cases[] = {
+        {cplx( 0,  0), cplx( 0,  0)},
+        {cplx( 1,  0), cplx( 1,  0)},
+        {cplx( 2,  0), cplx( 4,  0)},
+        {cplx(-1,  0), cplx( 1,  0)},
+        {cplx( 0,  1), cplx(-1,  0)},
+        {cplx( 0,  2), cplx(-4,  0)},
+        {cplx( 0, -1), cplx(-1,  0)},
+        {cplx( 1,  1), cplx( 0,  2)},
+        {cplx( 1, -1), cplx( 0, -2)},
+        {cplx(-1, -1), cplx( 0,  2)},
+        {cplx(0.5, 0), cplx(0.25, 0)},
+    };
+
+    const unsigned num_cases = sizeof(cases) / sizeof(test_case);
+    for (unsigned i = 0; i < num_cases; ++i)
+    {
+        const test_case& test = cases[i];
+        const std::complex<T> actual = std::__sqr(test.value);
+        assert(std::abs(actual.real() - test.expected.real()) < tolerance);
+        assert(std::abs(actual.imag() - test.expected.imag()) < tolerance);
+    }
+
+    const cplx nan1 = std::__sqr(cplx(NAN, 0));
+    assert(std::isnan(nan1.real()));
+    assert(std::isnan(nan1.imag()));
+
+    const cplx nan2 = std::__sqr(cplx(0, NAN));
+    assert(std::isnan(nan2.real()));
+    assert(std::isnan(nan2.imag()));
+
+    const cplx nan3 = std::__sqr(cplx(NAN, NAN));
+    assert(std::isnan(nan3.real()));
+    assert(std::isnan(nan3.imag()));
+
+    const cplx inf1 = std::__sqr(cplx(INFINITY, 0));
+    assert(std::isinf(inf1.real()));
+    assert(inf1.real() > 0);
+
+    const cplx inf2 = std::__sqr(cplx(0, INFINITY));
+    assert(std::isinf(inf2.real()));
+    assert(inf2.real() < 0);
+}
+
+int main()
+{
+    test<float>();
+    test<double>();
+    test<long double>();
+}

test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp

for (unsigned i = 0; i < N; ++i)
         }
         else if (testcases[i].real() == 1 && testcases[i].imag() == 0)
         {
             assert(r.real() == 0);
             assert(!std::signbit(r.real()));
             assert(r.imag() == 0);
             assert(std::signbit(r.imag()) == std::signbit(testcases[i].imag()));
         }
+        else if (testcases[i].real() == -1 && testcases[i].imag() == 0)
+        {
+            assert(r.real() == 0);
+            assert(!std::signbit(r.real()));
+            if (std::signbit(testcases[i].imag()))
+                is_about(r.imag(), -pi);
+            else
+                is_about(r.imag(),  pi);
+        }
         else if (std::isfinite(testcases[i].real()) && std::isinf(testcases[i].imag()))
         {
             assert(std::isinf(r.real()));
             assert(r.real() > 0);
             if (std::signbit(testcases[i].imag()))
                 is_about(r.imag(), -pi/2);
             else
                 is_about(r.imag(),  pi/2);

test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp

void test_edges()
     for (unsigned i = 0; i < N; ++i)
     {
         std::complex<double> r = asinh(testcases[i]);
         if (testcases[i].real() == 0 && testcases[i].imag() == 0)
         {
             assert(std::signbit(r.real()) == std::signbit(testcases[i].real()));
             assert(std::signbit(r.imag()) == std::signbit(testcases[i].imag()));
         }
+        else if (testcases[i].real() == 0 && std::abs(testcases[i].imag()) == 1)
+        {
+            assert(r.real() == 0);
+            assert(std::signbit(testcases[i].imag()) == std::signbit(r.imag()));
+            if (std::signbit(testcases[i].imag()))
+                is_about(r.imag(), -pi/2);
+            else
+                is_about(r.imag(),  pi/2);
+        }
         else if (std::isfinite(testcases[i].real()) && std::isinf(testcases[i].imag()))
         {
             assert(std::isinf(r.real()));
             assert(std::signbit(testcases[i].real()) == std::signbit(r.real()));
             if (std::signbit(testcases[i].imag()))
                 is_about(r.imag(), -pi/2);
             else
                 is_about(r.imag(),  pi/2);