[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Mikhail Maltsev via Phabricator via cfe-commits]

This revision was automatically updated to reflect the committed changes.
Closed by commit rCXX325510: [libcxx] Improve accuracy of complex asinh and 
acosh (authored by miyuki, committed by ).

Repository:
  rCXX libc++

https://reviews.llvm.org/D41629

Files:
  include/complex
  test/libcxx/numerics/complex.number/__sqr.pass.cpp
  test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp
  test/std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp

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
+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
@@ -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: test/libcxx/numerics/complex.number/__sqr.pass.cpp
===
--- test/libcxx/numerics/complex.number/__sqr.pass.cpp
+++ test/libcxx/numerics/complex.number/__sqr.pass.cpp
@@ -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.
+//
+//===--===//
+
+// 
+
+// template
+//   complex
+//   __sqr(const complex& x);
+
+#include 
+#include 
+
+template 
+void
+test()
+{
+const T tolerance = std::is_same::value 

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Mikhail Maltsev via Phabricator via cfe-commits]

miyuki added a comment.

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


https://reviews.llvm.org/D41629



___
cfe-commits mailing list
cfe-commits@lists.llvm.org
http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Marshall Clow via Phabricator via cfe-commits]

mclow.lists accepted this revision.
mclow.lists added a comment.
This revision is now accepted and ready to land.

LGTM. Do you need me to commit this?


https://reviews.llvm.org/D41629



___
cfe-commits mailing list
cfe-commits@lists.llvm.org
http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Mikhail Maltsev via Phabricator via cfe-commits]

miyuki added a comment.
Herald added a subscriber: christof.

ping


https://reviews.llvm.org/D41629



___
cfe-commits mailing list
cfe-commits@lists.llvm.org
http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Mikhail Maltsev via Phabricator via cfe-commits]

miyuki updated this revision to Diff 133408.
miyuki added a comment.

Added a test for __sqr


https://reviews.llvm.org/D41629

Files:
  include/complex
  test/libcxx/numerics/complex.number/__sqr.pass.cpp
  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: test/libcxx/numerics/complex.number/__sqr.pass.cpp
===
--- /dev/null
+++ test/libcxx/numerics/complex.number/__sqr.pass.cpp
@@ -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.
+//
+//===--===//
+
+// 
+
+// template
+//   complex
+//   __sqr(const complex& x);
+
+#include 
+#include 
+
+template 
+void
+test()
+{
+const T tolerance = std::is_same::value ? 1.e-6 : 1.e-14;
+
+typedef std::complex 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 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();
+test();
+test();
+}
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

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Marshall Clow via Phabricator via cfe-commits]

mclow.lists added a comment.

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"


https://reviews.llvm.org/D41629



___
cfe-commits mailing list
cfe-commits@lists.llvm.org
http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Mikhail Maltsev via Phabricator via cfe-commits]

miyuki added a comment.

ping^3


https://reviews.llvm.org/D41629



___
cfe-commits mailing list
cfe-commits@lists.llvm.org
http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Mikhail Maltsev via Phabricator via cfe-commits]

miyuki added a comment.

ping^2


https://reviews.llvm.org/D41629



___
cfe-commits mailing list
cfe-commits@lists.llvm.org
http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Mikhail Maltsev via Phabricator via cfe-commits]

miyuki added a comment.

ping


https://reviews.llvm.org/D41629



___
cfe-commits mailing list
cfe-commits@lists.llvm.org
http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Mikhail Maltsev via Phabricator via cfe-commits]

miyuki added inline comments.



Comment at: 
test/std/numerics/complex.number/complex.transcendentals/acosh.pass.cpp:59
+{
+assert(r.real() == 0);
+assert(!std::signbit(r.real()));

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.


https://reviews.llvm.org/D41629



___
cfe-commits mailing list
cfe-commits@lists.llvm.org
http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits

[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh

[Mikhail Maltsev via Phabricator via cfe-commits]

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
+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
@@ -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()) ==