This program demonstrates a bug in the complex sqrt function in gfortran. When complex number x = (0.0,-1.0), the result comes back as the negative of the expected result. Of course, the negative of the expected result is still a square root of the original number, but the f95 standard says that the principal square root should be returned, i.e. the real part of the result should be positive. complex x, y x = cmplx(0.0,-1.0) y = sqrt(x) write (*,*) x, y x = cmplx(0.0+1.0e-38,-1.0) y = sqrt(x) write (*,*) x, y x = cmplx(0.0-1.0e-38,-1.0) y = sqrt(x) write (*,*) x, y end The results I get are % gfortran -v Using built-in specs. Target: x86_64-unknown-linux-gnu Configured with: ../gcc/configure --prefix=/var/tmp/gfortran-20051007/irun --enable-languages=c,f95 Thread model: posix gcc version 4.1.0 20051007 (experimental) % gfortran sqrt.f % ./a.out ( 0.000000 , -1.000000 ) (-0.7071068 , 0.7071068 ) ( 9.9999994E-39, -1.000000 ) ( 0.7071068 ,-0.7071068 ) (-9.9999994E-39, -1.000000 ) ( 0.7071068 ,-0.7071068 ) Using a tiny number in place of the zero real part shows the correct results. This also affects the double complex version of sqrt.
Hmm, I think this is a bug in glibc's csqrt. as I get the following on powerpc-darwin: ( 0.000000 , -1.000000 ) ( 0.7071068 ,-0.7071068 ) ( 9.9999994E-39, -1.000000 ) ( 0.7071068 ,-0.7071068 ) (-9.9999994E-39, -1.000000 ) ( 0.7071068 ,-0.7071068 )
This is a bug in glibc's csqrt as gfortran just calls it. Please report it to them.
I don't think this is a glibc problem. gfortran will use a GCC builtin function. See gcc/fortran/mathbuiltins.def To confirm this, I get troutmask:kargl[210] ./z ( 0.000000 , -1.000000 ) (-0.7071068 , 0.7071068 ) ( 1.1754944E-38, -1.000000 ) ( 0.7071068 ,-0.7071068 ) (-1.1754944E-38, -1.000000 ) ( 0.7071068 ,-0.7071068 ) on amd64-*-freebsd, which does not have a csqrtf in libm. Note, I've replaced Mick's subnormal number with TINY(1.).
This is definitely a gcc bug. gfortran -fdump-tree-original csqrt.f yields (with removal of IO crap and shorting of long constants) MAIN__ () { complex4 y; complex4 x; x = __complex__ (0.0, -1.0e+0); y = __builtin_csqrtf (x); x = __complex__ (1.17549435082228750796873653722224e-38, -1.0e+0); y = __builtin_csqrtf (x); x = __complex__ (-1.175494350822287507968736537222e-38, -1.0e+0); y = __builtin_csqrtf (x); } Note the use of __builtin_csqrtf().
(In reply to comment #4) > Note the use of __builtin_csqrtf(). But that just calls into csqrt which is part of libm. Look at the asm output.
FreeBSD *does not* have a csqrt or csqrtf in libm! The problem is in intrinsics/c99_functions.c. I'm testing a patch.
(In reply to comment #6) > FreeBSD *does not* have a csqrt or csqrtf in libm! The problem is > in intrinsics/c99_functions.c. I'm testing a patch. But if csqrtf in c99_functions.c was pulled from glibc, then there is still a bug in glibc as this was reported with the target being x86_64-linux-gnu.
Yes, I agree glibc has the bug. But, libgfortran also has the bug and will be present on every OS that does not have a libm with csqrt and csqrtf. I have a patch for libgfortran.
Patch for libgfortran http://gcc.gnu.org/ml/gcc-patches/2005-10/msg00626.html
Subject: Bug 24313 CVSROOT: /cvs/gcc Module name: gcc Changes by: kargl@gcc.gnu.org 2005-10-11 23:35:27 Modified files: libgfortran : ChangeLog libgfortran/intrinsics: c99_functions.c gcc/testsuite : ChangeLog Added files: gcc/testsuite/gfortran.dg: csqrt_2.f Log message: PR libgfortran/24313 * c99_functions.c (csqrtf, csqrt): Fix choice of branch cut. Note csqrt{f} were imported from glibc, and this bug is still present there. glibc PR is 1146. Patches: http://gcc.gnu.org/cgi-bin/cvsweb.cgi/gcc/libgfortran/ChangeLog.diff?cvsroot=gcc&r1=1.320&r2=1.321 http://gcc.gnu.org/cgi-bin/cvsweb.cgi/gcc/libgfortran/intrinsics/c99_functions.c.diff?cvsroot=gcc&r1=1.16&r2=1.17 http://gcc.gnu.org/cgi-bin/cvsweb.cgi/gcc/gcc/testsuite/gfortran.dg/csqrt_2.f.diff?cvsroot=gcc&r1=NONE&r2=1.1 http://gcc.gnu.org/cgi-bin/cvsweb.cgi/gcc/gcc/testsuite/ChangeLog.diff?cvsroot=gcc&r1=1.6170&r2=1.6171
*** Bug 25017 has been marked as a duplicate of this bug. ***
(In reply to comment #9) > Patch for libgfortran > > http://gcc.gnu.org/ml/gcc-patches/2005-10/msg00626.html > It is a useful strategy to use libm but... It seems that glibc has still problems with csqrt. See also http://sources.redhat.com/bugzilla/show_bug.cgi?id=2181 http://sources.redhat.com/bugzilla/show_bug.cgi?id=2182 . Therefore it may be useful to check in intrinsics/configure if libm is broken in the architecture used, and if so use the code in intrinsics/c99_functions.c . A possible solution for csqrt is shown in http://www-zeuthen.desy.de/~hvogt/gfortran/acinclude.m4.diff . It is better to be independend from glibc coming with the OS. This enables to install gcc4(gfortran) without patching libm if it has bugs.
Subject: Bug 24313 Author: burnus Date: Sun Oct 29 11:56:56 2006 New Revision: 118142 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=118142 Log: 2006-10-28 Tobias Burnus <burnus@net-b.de> PR libgfortran/24313 * gfortran.dg/csqrt_2.f: Remove xfail *-*-linux-gnu. Modified: trunk/gcc/testsuite/ChangeLog trunk/gcc/testsuite/gfortran.dg/csqrt_2.f
For completeness, the two problems are in the meanwhile fixed in glibc: - csqrt bug, fixed 2005-10-13, http://sources.redhat.com/bugzilla/show_bug.cgi?id=1466 - cacosh bug, fixed 2006-08-03, http://sources.redhat.com/bugzilla/show_bug.cgi?id=2182