maxjhandsome
/
llvm-project
forked from OSchip/llvm-project
1
0
Fork 0
Code Issues Pull Requests Packages Releases Wiki Activity

[libcxx] Improve accuracy of complex asinh and acosh 

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.

Reviewers: EricWF, mclow.lists

Reviewed By: mclow.lists

Subscribers: christof, cfe-commits

Differential Revision: https://reviews.llvm.org/D41629

llvm-svn: 325510
This commit is contained in:
Mikhail Maltsev 2018-02-19 15:41:36 +00:00
parent 9b2aa42f00
commit 477b5f688c
4 changed files with 113 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
libcxx/include/complex
View File

@ -1125,6 +1125,17 @@ pow(const _Tp& __x, const complex<_Up>& __y)
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>
@ -1150,7 +1161,7 @@ asinh(const complex<_Tp>& __x)
}
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()));
}
@ -1184,7 +1195,7 @@ acosh(const complex<_Tp>& __x)
}
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()));
}
@ -1318,7 +1329,7 @@ acos(const complex<_Tp>& __x)
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()));

81
libcxx/test/libcxx/numerics/complex.number/__sqr.pass.cpp Normal file
View File

@ -0,0 +1,81 @@
//===----------------------------------------------------------------------===//
//
// 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>();
}

9
libcxx/test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp
View File

@ -54,6 +54,15 @@ void test_edges()
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()));

9
libcxx/test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp
View File

@ -44,6 +44,15 @@ void test_edges()
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()));