According to the C++ standard, acosh(complex<T>) should behave just like C99's cacosh(T complex) function. There, the branch cut is "at values less than 1 along the real axis" and the "range of a half-strip of non-negative values along the real axis and in the interval [-i*pi,+i*pi] along the imaginary axis." The implementation in tr1/complex gets this wrong. The result returned by __complex_acosh() are all wrong in the lower complex plain. It can be easily fixed. I'm attaching a patch.
Created attachment 25623 [details] patch to fix the bug
Thus, to understand and clarify why this has not been noticed so far, you are on a target which doesn't support in the underlying C library these complex functions, right? Because normally (eg, on Linux) these days we just forward to __builtin_cacosh*, the code you are touching is just a "surrogate", a "fallback", which doesn't get right all the special cases, NaNs, infinity. Anyway, a similar tweak would touch also the C++11 version in std:: Gaby, can you have a look to this, double check the patch? For your convenience the surrounding code is: template<typename _Tp> std::complex<_Tp> __complex_acosh(const std::complex<_Tp>& __z) { std::complex<_Tp> __t((__z.real() - __z.imag()) * (__z.real() + __z.imag()) - _Tp(1.0), _Tp(2.0) * __z.real() * __z.imag()); __t = std::sqrt(__t); return std::log(__t + __z); }
Ok, checked. Will commit the fix momentarily (4.6.3 too when the branch reopens)
Author: paolo Date: Thu Oct 27 11:00:25 2011 New Revision: 180563 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=180563 Log: 2011-10-27 Richard B. Kreckel <kreckel@ginac.de> Paolo Carlini <paolo.carlini@oracle.com> PR libstdc++/50880 * include/std/complex (__complex_acosh): Fix for __z.real() < 0. * include/tr1/complex (__complex_acosh): Likewise. * testsuite/26_numerics/complex/50880.cc: New. * testsuite/tr1/8_c_compatibility/complex/50880.cc: Likewise. Added: trunk/libstdc++-v3/testsuite/26_numerics/complex/50880.cc trunk/libstdc++-v3/testsuite/tr1/8_c_compatibility/complex/50880.cc Modified: trunk/libstdc++-v3/ChangeLog trunk/libstdc++-v3/include/std/complex trunk/libstdc++-v3/include/tr1/complex
On 10/27/2011 11:24 AM, paolo.carlini at oracle dot com wrote: > Thus, to understand and clarify why this has not been noticed so far, you are > on a target which doesn't support in the underlying C library these complex > functions, right? Because normally (eg, on Linux) these days we just forward to > __builtin_cacosh*, the code you are touching is just a "surrogate", a > "fallback", which doesn't get right all the special cases, NaNs, infinity. Well, I didn't "notice" it. Searching for bugs involving branch cut positions of inverse trig functions in a range of FOSS projects is a pet project of mine. ;-) BTW: I noticed that my patch tests if (__z.real() < -0.0), which is correct, but the sign of zero looks eccentric and might potentially confuse future readers. I suggest removing it. It doesn't matter at all, anyhow.
Indeed you are right about the sign, in terms at least of consistency with the rest of the fallback implementations which already have got quite a number of comparisons with zero with no special attention to its signedness (like '< _Tp()' or '> _Tp()'). I had already noticed that. As soon as I find a bit of time, we can also *consistently over all those cases* use __builtin_signbit, as suggested by Gaby elsewhere. I have to double check with the middle-end people that it doesn't generate library calls for the patch to be neat.
> As soon as I find a bit of > time, we can also *consistently over all those cases* use __builtin_signbit, as > suggested by Gaby elsewhere. I have to double check with the middle-end people > that it doesn't generate library calls for the patch to be neat. We also better double check whether the results stay correct. Thinking of it... Big Ooops! It turns out the patch makes it much better but still not entirely correct. On systems that feature a signed zero, it gives incorrect results when either * __z.real() is -0.0 and signbit(__z.imag()) * __z.real() < -1.0 and __z.imag() is -0.0 The first problem can be fixed by using signbit instead of -0.0, as you proposed, but the second needs more correction. The patch BC1 I'm attaching should fix these cases, too. But this is starting to look cumbersome. We are trying to construct a formula that is continuous except on the branch cut defined in C99. Why not just use the formula mentioned in CLTL2 [0] like in patch BC2 that I'm attaching? (The branch cuts of inverse trig functinos in C99 and Common Lisp are the same.) [0] <http://www-prod-gif.supelec.fr/docs/cltl/clm/node129.html#SECTION001653000000000000000>
Created attachment 25653 [details] BC1
Created attachment 25654 [details] BC2
Richard, I have no problems with BC2. This is code I wrote rather quickly a few years ago, adapting it from glibc, essentially, and then each year that went by, fewer and fewer systems used and tested it, because underlying C99 libc support is becoming often available (eg, for sure Linux and Darwin). Thus, please sleep on this, let's wait for comments from other people, and say, in a week or so we'll finalize the code for 4.7. Thanks again!
Paolo, I still intend to come forward with a patch for all these cases. Unfortunately, I was distracted by what I've just filed as Bug 50957.
In my opinion BC2 is fine, I can take of applying it, if you still endorse it.
(In reply to comment #2) > Thus, to understand and clarify why this has not been noticed so far, you are > on a target which doesn't support in the underlying C library these complex > functions, right? Because normally (eg, on Linux) these days we just forward to > __builtin_cacosh*, the code you are touching is just a "surrogate", a > "fallback", which doesn't get right all the special cases, NaNs, infinity. > > Anyway, a similar tweak would touch also the C++11 version in std:: > > Gaby, can you have a look to this, double check the patch? For your convenience > the surrounding code is: > > template<typename _Tp> > std::complex<_Tp> > __complex_acosh(const std::complex<_Tp>& __z) > { > std::complex<_Tp> __t((__z.real() - __z.imag()) > * (__z.real() + __z.imag()) - _Tp(1.0), > _Tp(2.0) * __z.real() * __z.imag()); > __t = std::sqrt(__t); > > return std::log(__t + __z); > } As I observed elsewhere, the test should be on the sign, no comparison against 0.0, so that signed zero is handled correctly.
(In reply to comment #9) > Created attachment 25654 [details] > BC2 Since we are talking about branch cut and prespectiving since zeros, I think we should avoid the the arithmetic z -/+ one, whee one is of a complex<T>. Rather the computation should be be directly on the components. This is to prevent signed zeros to have their mutated.
Ok, thanks for your feedback Gaby. Indeed, I also wondered if we shouldn't work with the components. Richard, can you send a version of Kahan's algorithm rewritten in terms of real and imag, I think we can all agree on that and resolve the PR for good.
Well, I guess this would be most of it: template<typename _Tp> std::complex<_Tp> __complex_acosh(const std::complex<_Tp>& __z) { return _Tp(2.0) * std::log(std::sqrt(_Tp(0.5) * (__z + _Tp(1.0))) + std::sqrt(_Tp(0.5) * (__z - _Tp(1.0)))); }
(In reply to comment #16) > Well, I guess this would be most of it: > > template<typename _Tp> > std::complex<_Tp> > __complex_acosh(const std::complex<_Tp>& __z) > { > return _Tp(2.0) * std::log(std::sqrt(_Tp(0.5) * (__z + _Tp(1.0))) > + std::sqrt(_Tp(0.5) * (__z - _Tp(1.0)))); > } looks good -- hoping for log implementation to do the right thing.
Author: paolo Date: Wed Nov 2 18:43:26 2011 New Revision: 180787 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=180787 Log: 2011-11-02 Richard B. Kreckel <kreckel@ginac.de> Paolo Carlini <paolo.carlini@oracle.com> PR libstdc++/50880 * include/std/complex (__complex_acosh): Fix in a better way, use Kahan's formula. * include/tr1/complex (__complex_acosh): Likewise. Modified: trunk/libstdc++-v3/include/std/complex trunk/libstdc++-v3/include/tr1/complex
Author: paolo Date: Wed Nov 2 18:43:42 2011 New Revision: 180788 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=180788 Log: 2011-11-02 Richard B. Kreckel <kreckel@ginac.de> Paolo Carlini <paolo.carlini@oracle.com> PR libstdc++/50880 * include/std/complex (__complex_acosh): Fix in a better way, use Kahan's formula. * include/tr1/complex (__complex_acosh): Likewise. Modified: trunk/libstdc++-v3/ChangeLog
Author: paolo Date: Wed Nov 2 21:54:24 2011 New Revision: 180804 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=180804 Log: 2011-11-02 Richard B. Kreckel <kreckel@ginac.de> Paolo Carlini <paolo.carlini@oracle.com> PR libstdc++/50880 * include/std/complex (__complex_acosh): Fix in a better way, use Kahan's formula. * include/tr1/complex (__complex_acosh): Likewise. Added: branches/gcc-4_6-branch/libstdc++-v3/testsuite/26_numerics/complex/50880.cc branches/gcc-4_6-branch/libstdc++-v3/testsuite/tr1/8_c_compatibility/complex/50880.cc Modified: branches/gcc-4_6-branch/libstdc++-v3/ChangeLog branches/gcc-4_6-branch/libstdc++-v3/include/std/complex branches/gcc-4_6-branch/libstdc++-v3/include/tr1/complex
Fixed for 4.6.3 and mainline.
(In reply to comment #16) > Well, I guess this would be most of it: > > template<typename _Tp> > std::complex<_Tp> > __complex_acosh(const std::complex<_Tp>& __z) > { > return _Tp(2.0) * std::log(std::sqrt(_Tp(0.5) * (__z + _Tp(1.0))) > + std::sqrt(_Tp(0.5) * (__z - _Tp(1.0)))); > } [Sorry for my temporary absence.] For future reference, some final remarks: * Yes, that is a good implementation for this "fallback". It relies on __z - _Tp(1.0) not "mutating" the sign of __z's imag part. * If the complex log does not do the right thing, all is lost anyways (besides, __complex_asinh relies on it, too). * My patch BC1 was flawed. It contains code trying to work around a ctor doing multiplication (fixed in PR48760) * My patch BC2 was flawed for the same reason: __z - __one doesn't preserve the sign of __z's imag part. Looks good. Thanks!
Thanks Richard for double checking!
By the way, if isn't clear already, I would be *really* curious to know which specific targets by now can't just enable the builtins, eg, their libc doesn't provide the C99 complex functions.
(In reply to comment #25) > By the way, if isn't clear already, I would be *really* curious to know which > specific targets by now can't just enable the builtins, eg, their libc doesn't > provide the C99 complex functions. Sadly, the libc functions aren't always correct either. But that is another story: http://sourceware.org/bugzilla/show_bug.cgi?id=13305