[Bug libstdc++/89991] Complex numbers: Calculation of imaginary part is not correct
kargl at gcc dot gnu.org
gcc-bugzilla@gcc.gnu.org
Mon Apr 8 21:46:00 GMT 2019
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991
--- Comment #19 from kargl at gcc dot gnu.org ---
(In reply to Steve Kargl from comment #18)
> On Mon, Apr 08, 2019 at 08:03:36PM +0000, pinskia at gcc dot gnu.org wrote:
> > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89991
> >
> > --- Comment #17 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
> > >Doesn't this gets the wrong answer for __y = -0 (as -0 < 0 is false)?
> >
> > No, you missed this part:
> > // The branch cut is on the negative axis.
>
> No, I didn't miss that part.
>
> > So maybe the bug is inside FreeBSD and Window's libm. Glibc fixed the branch
> > cuts issues back in 2012 for csqrt but the other OS's did not change theirs.
>
> For the C++ code in comment, on x86_64-*-freebsd.
>
> % g++8 -o z a.cpp -lm && ./z
> z = (-1.84250315177824e-07,-0)
> pow(z,0.5) = (2.62836076598358e-20,-0.000429243887758258)
> sqrt(z) = (0,0.000429243887758258)
> sqrt(conj(z)) = (0,0.000429243887758258)
> conj(sqrt(z)) = (0,-0.000429243887758258)
>
> The last two lines are definitely wrong.
>
> troutmask:sgk[209] nm z | grep csqrt
> troutmask:sgk[210] nm z | grep sqrt
> 000000000040156b W _ZSt14__complex_sqrtIdESt7complexIT_ERKS2_
> 000000000040143d W _ZSt4sqrtIdESt7complexIT_ERKS2_
> U sqrt@@FBSD_1.0
>
> There is no reference to csqrt in the exectuable. If I change
> /usr/local/lib/gcc8/include/c++/complex to use copysign
> to account for __y = -0 like
>
> template<typename _Tp>
> complex<_Tp>
> __complex_sqrt(const complex<_Tp>& __z)
> {
> _Tp __x = __z.real();
> _Tp __y = __z.imag();
>
> if (__x == _Tp())
> {
> _Tp __t = sqrt(abs(__y) / 2);
> // return complex<_Tp>(__t, __y < _Tp() ? -__t : __t);
> return complex<_Tp>(__t, copysign(__t, __y));
> }
> else
> {
> _Tp __t = sqrt(2 * (std::abs(__z) + abs(__x)));
> _Tp __u = __t / 2;
> // return __x > _Tp()
> // ? complex<_Tp>(__u, __y / __t)
> // : complex<_Tp>(abs(__y) / __t, __y < _Tp() ? -__u : __u);
> return __x > _Tp()
> ? complex<_Tp>(__u, __y / __t)
> : complex<_Tp>(abs(__y) / __t, copysign(__u, __y));
> }
> }
>
>
> The C++ code in comment #10 gives
>
> g++8 -o z a.cpp -lm && ./z
> z = (-1.84250315177824e-07,-0)
> pow(z,0.5) = (2.62836076598358e-20,-0.000429243887758258)
> sqrt(z) = (0,-0.000429243887758258)
> sqrt(conj(z)) = (0,0.000429243887758258)
> conj(sqrt(z)) = (0,0.000429243887758258)
>
> The correct answer. QED.
BTW, if change /usr/local/lib/gcc8/include/c++/complex back to
no using copysign(), and instead change
#if _GLIBCXX_USE_C99_COMPLEX
inline __complex__ float
__complex_sqrt(__complex__ float __z) { return __builtin_csqrtf(__z); }
to
#if _GLIBCXX_USE_C99_COMPLEX || SOMETHING_UGLY
inline __complex__ float
__complex_sqrt(__complex__ float __z) { return __builtin_csqrtf(__z); }
and do
g++8 -DSOMETHING_UGLY -o z a.cpp -lm && ./z
z = (-1.84250315177824e-07,-0)
pow(z,0.5) = (2.62836076598358e-20,-0.000429243887758258)
sqrt(z) = (0,-0.000429243887758258)
sqrt(conj(z)) = (0,0.000429243887758258)
conj(sqrt(z)) = (0,0.000429243887758258)
I get the expected. So, if you're on a system that has
_GLIBCXX_USE_C99_COMPLEX, you won't see the bug.
It is likely that everywhere that a construct of the
form __y < _Tp() ? -__u : __u appear, it needs to use
copysign.
More information about the Gcc-bugs
mailing list