miyuki created this revision.
miyuki added reviewers: EricWF, mclow.lists.
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.
https://reviews.llvm.org/D41629
Files:
include/complex
test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp
test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp
Index: test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp
===================================================================
--- test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp
+++ test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp
@@ -44,6 +44,15 @@
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()));
Index: test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp
===================================================================
--- test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp
+++ test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp
@@ -54,6 +54,15 @@
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()));
Index: include/complex
===================================================================
--- include/complex
+++ include/complex
@@ -1125,6 +1125,17 @@
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 @@
}
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 @@
}
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 @@
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()));
Index: test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp
===================================================================
--- test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp
+++ test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp
@@ -44,6 +44,15 @@
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()));
Index: test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp
===================================================================
--- test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp
+++ test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp
@@ -54,6 +54,15 @@
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()));
Index: include/complex
===================================================================
--- include/complex
+++ include/complex
@@ -1125,6 +1125,17 @@
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 @@
}
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 @@
}
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 @@
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()));
_______________________________________________
cfe-commits mailing list
cfe-commits@lists.llvm.org
http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits
