[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 Harald Anlauf changed: What|Removed |Added CC||anlauf at gmx dot de --- Comment #1 from Harald Anlauf --- I can reproduce the exception only for -O0, but not for -O1. Note that the exponent is real (1e0), not an integer, and x**y is singular for x=0 and arbitrary y. Did you expect that the result should be well-defined?
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #2 from Vittorio Zecca --- I believe -O0 is the default optimization level, so it is important. Of course the code I sent you is a fragment from a much larger program, kind of c**x with c complex and x real. It is not possible to make x into an integer. When x is zero and y is real, x**y is singular for y<=0, right? Also, I do not understand why I get SIGFPE on either zero or invalid -ffpe-trap suboption, but this is a lesser issue.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #3 from Harald Anlauf --- (In reply to Vittorio Zecca from comment #2) > I believe -O0 is the default optimization level, so it is important. > > Of course the code I sent you is a fragment from a much larger program, > kind of c**x with c complex and x real. It is not possible to make x > into an integer. > > When x is zero and y is real, x**y is singular for y<=0, right? I think you are missing the intricacies of complex arithmetics. x**y = exp (y * log(x)) has an *essential singularity* at x=0, which is the starting point of a branch cut (usually the negative real axis). See also cpow(3). The man page for pow(3) covers only real arithmetics and does not apply. > Also, I do not understand why I get SIGFPE on either zero or invalid > -ffpe-trap suboption, > but this is a lesser issue. A quick search did not turn up any result whether IEEE specifies a defined behavior for your case. Maybe you are more successful. I also could not find anything in the Fortran standard.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #4 from Vittorio Zecca --- I am happy to refresh my complex analysis study of long ago. The singularity of log(x) in zero is not essential. Essential singularity means the Laurent seriesis infinite in both directions? z**-k and z**+k roughly, but log(z) Laurent series is infinite only for z**+k. I hope to remember correctly. But exp(y*log(x)) may well be essential, however. Anyway I am not sure this is an essential (forgive the pun) reason to raise an exception Also I do not understand the different behaviour with different levels of optimization, and I confirm the a.out executable runs fine under valgrind. And the code runs fine with Intel ifort. And I really do not see how complex zero raised to a positive power should raise an exception. 2013/6/29 anlauf at gmx dot de > http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 > > --- Comment #3 from Harald Anlauf --- > (In reply to Vittorio Zecca from comment #2) > > I believe -O0 is the default optimization level, so it is important. > > > > Of course the code I sent you is a fragment from a much larger program, > > kind of c**x with c complex and x real. It is not possible to make x > > into an integer. > > > > When x is zero and y is real, x**y is singular for y<=0, right? > > I think you are missing the intricacies of complex arithmetics. > > x**y = exp (y * log(x)) has an *essential singularity* at x=0, > which is the starting point of a branch cut (usually the negative > real axis). See also cpow(3). > > The man page for pow(3) covers only real arithmetics and does not apply. > > > Also, I do not understand why I get SIGFPE on either zero or invalid > > -ffpe-trap suboption, > > but this is a lesser issue. > > A quick search did not turn up any result whether IEEE specifies > a defined behavior for your case. Maybe you are more successful. > I also could not find anything in the Fortran standard. > > -- > You are receiving this mail because: > You reported the bug. >
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #5 from Harald Anlauf --- (In reply to Vittorio Zecca from comment #4) > I am happy to refresh my complex analysis study of long ago. > The singularity of log(x) in zero is not essential. Right. > Essential singularity means the Laurent seriesis infinite in both > directions? > z**-k and z**+k roughly, but log(z) Laurent series is infinite only for > z**+k. > I hope to remember correctly. > But exp(y*log(x)) may well be essential, however. Yes, since exp(z) has an essential singularity at complex infinity. > Anyway I am not sure this is an essential (forgive the pun) reason to raise > an exception So what should the correct result be? > Also I do not understand the different behaviour with different levels of > optimization, I think that compile-time optimization realizes that the exponent y is actually exactly a positive integer and does some simplification. At -O0, you get an evaluation by the run-time library. > and I confirm the a.out executable runs fine under valgrind. > And the code runs fine with Intel ifort. > And I really do not see how complex zero raised to a positive power should > raise an exception. Well, you actually provide a non-integer (real or complex) exponent, even if it is accidentally a positive integer.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749
--- Comment #6 from Harald Anlauf ---
(In reply to Harald Anlauf from comment #5)
> (In reply to Vittorio Zecca from comment #4)
> > Also I do not understand the different behaviour with different levels of
> > optimization,
>
> I think that compile-time optimization realizes that the exponent y
> is actually exactly a positive integer and does some simplification.
> At -O0, you get an evaluation by the run-time library.
Looking at the dump tree, one find that the following is generated:
D.1825 = __builtin_cpowf (D.1824, __complex__ (1.0e+0, 0.0));
Without optimization, this is converted into a call to cpowf.
A straigtforward C example shows that the exception is generated
within libm:
#define _GNU_SOURCE
#include
#include
#include
main()
{
complex x, y, z;
x = 0.;
y = 1.;
feenableexcept (FE_DIVBYZERO | FE_OVERFLOW | FE_INVALID);
z = cpowf (x, y);
printf("z = %f + %f * i\n", creal(z), cimag(z));
}
(gdb) run
Starting program: /home/anlauf/a.out
Program received signal SIGFPE, Arithmetic exception.
0xb7f98ea3 in clogf () from /lib/libm.so.6
Missing separate debuginfos, use: zypper install
glibc-debuginfo-2.14.1-14.27.1.i686
(gdb) where
#0 0xb7f98ea3 in clogf () from /lib/libm.so.6
#1 0xb7f9a4d4 in cpowf () from /lib/libm.so.6
#2 0x080485a2 in main () at check-cpow.c:12
That's with glibc-2.14.1.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #7 from Vittorio Zecca --- Looking at the source code for cpowf, it computes x**y straightly as exp(y*log(x)), no check on y or x. I am not happy with that, because I expect that complex zero raised to any (real or integer) positive power should be zero. No exceptions raised. I do not know if this is correct, but it looks reasonable to me. Computing complex zero**1 delivers (0.0,0.0), but complex zero**1.0 raises an exception. How strange. This is reading the base and exponent, no optimization involved. As in: complex zero read *,zero,i,x print *,zero,i,x print *,zero**i print *,zero**x When both numbers are complex, the standards (F95, F2003, and F2008) state that "the value of the operation x1**x2 is the principal value of x1 raised to x2". Does it help? When aimag(x2) is zero what is the principal value? I cannot believe that zero**1e0 is a singularity. The Intel ifort compiler delivers the "reasonable" result, it even goes beyond that and computes zero**zero as (1.0,0.0)! And I still believe the result should not depend on the optimization level. Note that compiling without -ffpe-trap the result with default optimization, -O0, is ( 0., -0.), higher optimizations deliver ( 0., 0.). Another mystery.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #8 from Dominique d'Humieres --- > And I still believe the result should not depend on the optimization level. Well, it does and IMO it has to (see below). > Note that compiling without -ffpe-trap the result > with default optimization, -O0, is ( 0., -0.), > higher optimizations deliver ( 0., 0.). > Another mystery. IMO there is no mystery. zero is a constant that is propagated during the optimization pass: looking at pr57749.f90.165t.optimized, I see D.1881 = __complex__ (0.0, 0.0); _gfortran_transfer_complex_write (&dt_parm.0, &D.1881, 4); i.e., zero**1e0 is computed during the optimization (hence cannot generate exception at run time). In pr57749.f90.003t.original, I see D.1881 = __builtin_cpowf (D.1880, __complex__ (1.0e+0, 0.0)); _gfortran_transfer_complex_write (&dt_parm.0, &D.1881, 4); so at run time you call a "buggy(?)" cpowf that generates the exception and (0.0,-0.0). As said in comment #6, this is not a gfortran bug, but a libgcc/libm one to be reported upstrean.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #9 from Vittorio Zecca --- Yes, I agree that there is a bug, and IMO it is in cpow/cpowf/cpowl. With -ffpe-trap=invalid,zero, as I wrote earlier, complex zero**I where I is integer equal to one, does not raise any exception and delivers the expected (by me) result, zero. Complex zero**X where X is real equal to one raises exception. Complex zero**C where C is complex equal to (1e0,0e0) raises exception. This is because glibc computes x**y as exp(y*log(x)) so it fails when x=0, but IMO it should not compute log(x) if x=0 and y>0 (y real), it should take a shortcut instead and deliver complex zero.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #10 from Dominique d'Humieres --- > Yes, I agree that there is a bug, ... Then you should report to the library maintainers, not to gfortran.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #11 from Harald Anlauf --- (In reply to Vittorio Zecca from comment #9) > Yes, I agree that there is a bug, and IMO it is in cpow/cpowf/cpowl. Vittorio, I'm even still not convinced that there is a bug here, you might be hitting undefined behavior. > With -ffpe-trap=invalid,zero, > as I wrote earlier, complex zero**I where I is integer equal to one, > does not raise > any exception and delivers the expected (by me) result, zero. > Complex zero**X where X is real equal to one raises exception. > Complex zero**C where C is complex equal to (1e0,0e0) raises exception. > This is because glibc computes x**y as exp(y*log(x)) so it fails when x=0, > but IMO it should not compute log(x) if x=0 and y>0 (y real), > it should take a shortcut instead and deliver complex zero. Now you're making several assumptions that you have to check. As I tried to explain, x ** y = exp (y * log(x)) is a rather complex (now there's also a pun intended) function: Consider that x and y are both complex. For x=0, it is *not* an analytic function of y at y=1. For arbitrary complex y, it is *not* an analytic function of x at x=0. Now you have to find out: - Is cpowf(x,y) an analytic function of both x and y? Then you're moving on thin ice and possibly invoke undefined behavior. - Is cpowf(x,y) identical to Fortran's x ** y for complex x, y? See above. - Does IEEE or the Fortran standard require what the result should be? If not, the operation is undefined. You are too much confined on assuming that y is real. Does the Fortran standard separate complex ** real from complex ** complex? I think not. Neither does C. I know that one define suitable paths in the complex plane that produce the result you desire for y real, but I am not convinced that this special case is covered by a C or Fortran standard. I'd recommend to fix your program to avoid numerical singularities of the above type (either change the exponent to integer, or avoid base zero). Or complain to the libm maintainers.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #12 from Vittorio Zecca --- > --- Comment #10 from Dominique d'Humieres --- >> Yes, I agree that there is a bug, ... > > Then you should report to the library maintainers, not to gfortran. The notion that this may be a library bug was not obvious to me from the beginning. I think that if the gfortran developers suspect there is a bug in a library called by gfortran to handle a numeric intrinsic operator like exponentiation, then they should be worried, more than I am.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #13 from Dominique d'Humieres --- We all whould have read the manual first: >From 'man cpow' ... DESCRIPTION cpow(x, y) returns the complex number x raised to the complex power y. cpow(x,y) is equivalent to cexp(y * clog(x)). As such, cpow(x, y) has a branch cut along the negative real axis for the first argument, and the equality cpow(conj(x),conj(y)) = conj(cpow(x,y)) holds for all x and y. SPECIAL VALUES For special values, see clog and cexp. ... >From 'man clog' ... SPECIAL VALUES The conjugate symmetry of clog() is used to abbreviate the specification of special values. clog(-0 + 0i) returns -inf + Pi i and raises the divide-by-zero flag. clog(0 + 0i) returns -inf + 0i and raises the divide-by-zero flag. ^^ ... >From 'man cexp' ... SPECIAL VALUES ... For the following two cases, cis(y) denotes cos(y) + I*sin(y). cexp(-inf + yi) returns 0*cis(y), for finite y. ... So this PR is not a bug, but a documented feature, and should be closed as INVALID.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #14 from Andreas Schwab --- C11 G.6.4.1 The cpow functions The cpow functions raise floating-point exceptions if appropriate for the calculation of the parts of the result, and may also raise spurious floating-point exceptions.379) 379) This allows cpow(z, c) to be implemented as cexp(c clog(z)) without precluding implementations that treat special cases more carefully.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 kargl at gcc dot gnu.org changed: What|Removed |Added Status|UNCONFIRMED |RESOLVED CC||kargl at gcc dot gnu.org Resolution|--- |INVALID --- Comment #15 from kargl at gcc dot gnu.org --- 1) This is not a gfortran issue as it has been shown to be a glibc issue. 2) It would be completely idiotic to pessimize the implementation of x**y by implementing some complex sequence of IF statements for a large number of special values. This is especially true for an OS that uses an alternate libm, which may already handle the special cases. 3) If one is writing portable Fortran code where this special case can occur, then it would be quite stupid to not test for this special case.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749
--- Comment #16 from Vittorio Zecca ---
You are being a little too hard on me, but so be it.
I believe there is only one special case, base==0,
and that there are only two ifs to put in cpow to avoid the floating exception
and give the expected result(I am simplifying here, also because I do
not use C):
if(base==0)
{
if(exponent>0) return 0; else raise hell;
}
The actual code where the original issue occurred had the exponentiation
in the deep of nested loops, it would have been rather time consuming
to test base==0
at the Fortran level
And I still do not understand why if the exponent is integer no
exception is raised and
the expected result zero is delivered.
As in the following fragment (with option -ffpe-trap=zero,invalid):
complex x
x=cmplx(0e0,0e0)
i=2
r=2e0
print *,x**i ! no exception raised delivers zero
print *,x**r ! exception raised
end
The Intel ifort and NAG nagfor compilers raise no exceptions and
deliver the expected result.
With my best wishes of good work to everybody.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #17 from Dominique d'Humieres --- > The actual code where the original issue occurred had the exponentiation > in the deep of nested loops, it would have been rather time consuming > to test base==0 > at the Fortran level Sorry, but I don't buy the argument. The cost of the test will likely be the same in gfortran or in the library. > And I still do not understand why if the exponent is integer no > exception is raised and > the expected result zero is delivered. When the exponent is integer, the computation is done through multiplications and zero times anything finite is zero.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749
--- Comment #18 from Steve Kargl ---
On Tue, Jul 02, 2013 at 07:21:54AM +, zeccav at gmail dot com wrote:
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749
>
> --- Comment #16 from Vittorio Zecca ---
> and that there are only two ifs to put in cpow to avoid the floating exception
> and give the expected result(I am simplifying here, also because I do
> not use C):
>
> if(base==0)
> {
> if(exponent>0) return 0; else raise hell;
> }
>
> The actual code where the original issue occurred had the exponentiation
> in the deep of nested loops, it would have been rather time consuming
> to test base==0
> at the Fortran level
It will be time consuming wherever it is tested. That's my
entire point and why the C11 standard permits cpow(z,c) to
be implemented as cexp(c*clog(z)) without trying to deal
with all of the special cases (or accuracy issues).
> And I still do not understand why if the exponent is integer no
> exception is raised and
> the expected result zero is delivered.
> As in the following fragment (with option -ffpe-trap=zero,invalid):
> complex x
> x=cmplx(0e0,0e0)
> i=2
> r=2e0
> print *,x**i ! no exception raised delivers zero
The compiler knows that i is an integer, and the above case
it is 2. The compiler evaluates x**2 as x*x.
> print *,x**r ! exception raised
> end
> The Intel ifort and NAG nagfor compilers raise no exceptions and
> deliver the expected result.
While it may be possible for a compiler to determine that r in
r=2e0 is an integral value of 2 and replace x**r by x*x, I suspect
that it will fail in the general case. What does
r = 8.125
x = (0.,0.)
print *, x**r ! x**8 * sqrt(sqrt(sqrt(x))) = 0.
or
r = 1. / 3
x = (0.,0,)
print *, x**r ! cube root of x?
produce?
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749
--- Comment #19 from Harald Anlauf ---
(In reply to Vittorio Zecca from comment #16)
> You are being a little too hard on me, but so be it.
>
> I believe there is only one special case, base==0,
> and that there are only two ifs to put in cpow to avoid the floating
> exception
> and give the expected result(I am simplifying here, also because I do
> not use C):
>
> if(base==0)
> {
> if(exponent>0) return 0; else raise hell;
> }
>
> The actual code where the original issue occurred had the exponentiation
> in the deep of nested loops, it would have been rather time consuming
> to test base==0
> at the Fortran level
>
>
> And I still do not understand why if the exponent is integer no
> exception is raised and
> the expected result zero is delivered.
> As in the following fragment (with option -ffpe-trap=zero,invalid):
> complex x
> x=cmplx(0e0,0e0)
> i=2
> r=2e0
> print *,x**i ! no exception raised delivers zero
> print *,x**r ! exception raised
> end
> The Intel ifort and NAG nagfor compilers raise no exceptions and
> deliver the expected result.
> With my best wishes of good work to everybody.
You obviously haven't tried other compilers.
With xlf, the result also depends on compiler flags:
Either (0.,0.), (NaNQ,NaNQ), or
Trace/BPT trap (core dumped)
I think you should accept that your code invokes undefined behavior
and needs fixing.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #20 from Vittorio Zecca --- I did try Intel ifort and NAG nagfor, as I already wrote, and also Solaris sunf90, and all of them accept complex zero raised to an integer or real power greater than zero, without raising any exception end delivering zero, the result I expected. It also works with a complex exponent whose real part is greater than zero, and zero imaginary part. I have not xlf, so I could not try it. Also, I understand the responsibility for the exponentiation behaviour is in glibc, not gfortran, so we may well close the discussion here.
[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749 --- Comment #21 from Vittorio Zecca --- If the cost is the same, this is a good reason to let the library handle this special case and let the programmer think about his/her business. Shouldn't system software ease the life of programmers? Shouldn't gfortran be consistent and raise the exception all the time, with integer real and complex exponent, or never? I mean, if the complex exponentiation has a (essential) singularity there, this is regardless of the Fortran type of the exponent. About portability, on my AMD 64 bit CPU with Fedora 18 I confirm that Intel ifort, NAG nagfor, and Solaris sunf95 do not raise any exception and deliver zero. A Fortran programmer may not know C and may think that a bad ** is a Fortran implementation issue. Anyway, I am preparing a note for the glibc people.
