Skip to content

Commit 477b5f6

Browse files
miyukimiyuki
committedFeb 19, 2018
[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
·
llvmorg-21-initllvmorg-7.0.0-rc1
1 parent 9b2aa42 commit 477b5f6

File tree

4 files changed

+113
-3
lines changed

4 files changed

+113
-3
lines changed
 

‎libcxx/include/complex

+14-3
Original file line numberDiff line numberDiff line change
@@ -1125,6 +1125,17 @@ pow(const _Tp& __x, const complex<_Up>& __y)
11251125
return _VSTD::pow(result_type(__x), result_type(__y));
11261126
}
11271127

1128+
// __sqr, computes pow(x, 2)
1129+
1130+
template<class _Tp>
1131+
inline _LIBCPP_INLINE_VISIBILITY
1132+
complex<_Tp>
1133+
__sqr(const complex<_Tp>& __x)
1134+
{
1135+
return complex<_Tp>((__x.real() - __x.imag()) * (__x.real() + __x.imag()),
1136+
_Tp(2) * __x.real() * __x.imag());
1137+
}
1138+
11281139
// asinh
11291140

11301141
template<class _Tp>
@@ -1150,7 +1161,7 @@ asinh(const complex<_Tp>& __x)
11501161
}
11511162
if (__libcpp_isinf_or_builtin(__x.imag()))
11521163
return complex<_Tp>(copysign(__x.imag(), __x.real()), copysign(__pi/_Tp(2), __x.imag()));
1153-
complex<_Tp> __z = log(__x + sqrt(pow(__x, _Tp(2)) + _Tp(1)));
1164+
complex<_Tp> __z = log(__x + sqrt(__sqr(__x) + _Tp(1)));
11541165
return complex<_Tp>(copysign(__z.real(), __x.real()), copysign(__z.imag(), __x.imag()));
11551166
}
11561167

@@ -1184,7 +1195,7 @@ acosh(const complex<_Tp>& __x)
11841195
}
11851196
if (__libcpp_isinf_or_builtin(__x.imag()))
11861197
return complex<_Tp>(abs(__x.imag()), copysign(__pi/_Tp(2), __x.imag()));
1187-
complex<_Tp> __z = log(__x + sqrt(pow(__x, _Tp(2)) - _Tp(1)));
1198+
complex<_Tp> __z = log(__x + sqrt(__sqr(__x) - _Tp(1)));
11881199
return complex<_Tp>(copysign(__z.real(), _Tp(0)), copysign(__z.imag(), __x.imag()));
11891200
}
11901201

@@ -1318,7 +1329,7 @@ acos(const complex<_Tp>& __x)
13181329
return complex<_Tp>(__pi/_Tp(2), -__x.imag());
13191330
if (__x.real() == 0 && (__x.imag() == 0 || isnan(__x.imag())))
13201331
return complex<_Tp>(__pi/_Tp(2), -__x.imag());
1321-
complex<_Tp> __z = log(__x + sqrt(pow(__x, _Tp(2)) - _Tp(1)));
1332+
complex<_Tp> __z = log(__x + sqrt(__sqr(__x) - _Tp(1)));
13221333
if (signbit(__x.imag()))
13231334
return complex<_Tp>(abs(__z.imag()), abs(__z.real()));
13241335
return complex<_Tp>(abs(__z.imag()), -abs(__z.real()));

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

+81
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
//===----------------------------------------------------------------------===//
2+
//
3+
// The LLVM Compiler Infrastructure
4+
//
5+
// This file is dual licensed under the MIT and the University of Illinois Open
6+
// Source Licenses. See LICENSE.TXT for details.
7+
//
8+
//===----------------------------------------------------------------------===//
9+
10+
// <complex>
11+
12+
// template<class T>
13+
// complex<T>
14+
// __sqr(const complex<T>& x);
15+
16+
#include <complex>
17+
#include <cassert>
18+
19+
template <class T>
20+
void
21+
test()
22+
{
23+
const T tolerance = std::is_same<T, float>::value ? 1.e-6 : 1.e-14;
24+
25+
typedef std::complex<T> cplx;
26+
struct test_case
27+
{
28+
cplx value;
29+
cplx expected;
30+
};
31+
32+
const test_case cases[] = {
33+
{cplx( 0, 0), cplx( 0, 0)},
34+
{cplx( 1, 0), cplx( 1, 0)},
35+
{cplx( 2, 0), cplx( 4, 0)},
36+
{cplx(-1, 0), cplx( 1, 0)},
37+
{cplx( 0, 1), cplx(-1, 0)},
38+
{cplx( 0, 2), cplx(-4, 0)},
39+
{cplx( 0, -1), cplx(-1, 0)},
40+
{cplx( 1, 1), cplx( 0, 2)},
41+
{cplx( 1, -1), cplx( 0, -2)},
42+
{cplx(-1, -1), cplx( 0, 2)},
43+
{cplx(0.5, 0), cplx(0.25, 0)},
44+
};
45+
46+
const unsigned num_cases = sizeof(cases) / sizeof(test_case);
47+
for (unsigned i = 0; i < num_cases; ++i)
48+
{
49+
const test_case& test = cases[i];
50+
const std::complex<T> actual = std::__sqr(test.value);
51+
assert(std::abs(actual.real() - test.expected.real()) < tolerance);
52+
assert(std::abs(actual.imag() - test.expected.imag()) < tolerance);
53+
}
54+
55+
const cplx nan1 = std::__sqr(cplx(NAN, 0));
56+
assert(std::isnan(nan1.real()));
57+
assert(std::isnan(nan1.imag()));
58+
59+
const cplx nan2 = std::__sqr(cplx(0, NAN));
60+
assert(std::isnan(nan2.real()));
61+
assert(std::isnan(nan2.imag()));
62+
63+
const cplx nan3 = std::__sqr(cplx(NAN, NAN));
64+
assert(std::isnan(nan3.real()));
65+
assert(std::isnan(nan3.imag()));
66+
67+
const cplx inf1 = std::__sqr(cplx(INFINITY, 0));
68+
assert(std::isinf(inf1.real()));
69+
assert(inf1.real() > 0);
70+
71+
const cplx inf2 = std::__sqr(cplx(0, INFINITY));
72+
assert(std::isinf(inf2.real()));
73+
assert(inf2.real() < 0);
74+
}
75+
76+
int main()
77+
{
78+
test<float>();
79+
test<double>();
80+
test<long double>();
81+
}

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

+9
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,15 @@ void test_edges()
5454
assert(r.imag() == 0);
5555
assert(std::signbit(r.imag()) == std::signbit(testcases[i].imag()));
5656
}
57+
else if (testcases[i].real() == -1 && testcases[i].imag() == 0)
58+
{
59+
assert(r.real() == 0);
60+
assert(!std::signbit(r.real()));
61+
if (std::signbit(testcases[i].imag()))
62+
is_about(r.imag(), -pi);
63+
else
64+
is_about(r.imag(), pi);
65+
}
5766
else if (std::isfinite(testcases[i].real()) && std::isinf(testcases[i].imag()))
5867
{
5968
assert(std::isinf(r.real()));

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

+9
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,15 @@ void test_edges()
4444
assert(std::signbit(r.real()) == std::signbit(testcases[i].real()));
4545
assert(std::signbit(r.imag()) == std::signbit(testcases[i].imag()));
4646
}
47+
else if (testcases[i].real() == 0 && std::abs(testcases[i].imag()) == 1)
48+
{
49+
assert(r.real() == 0);
50+
assert(std::signbit(testcases[i].imag()) == std::signbit(r.imag()));
51+
if (std::signbit(testcases[i].imag()))
52+
is_about(r.imag(), -pi/2);
53+
else
54+
is_about(r.imag(), pi/2);
55+
}
4756
else if (std::isfinite(testcases[i].real()) && std::isinf(testcases[i].imag()))
4857
{
4958
assert(std::isinf(r.real()));

0 commit comments

Comments
 (0)
Please sign in to comment.