[PATCH] D41629: [libcxx] Improve accuracy of complex asinh and acosh
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
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
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
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
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
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
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
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
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
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
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()) ==