[Bug fortran/57749] New: -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0

[Bug fortran/57749] New: -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0

[zeccav at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

            Bug ID: 57749
           Summary: -ffpe-trap=zero or invalid produces SIGFPE on complex
                    zero ** 1e0
           Product: gcc
           Version: 4.8.1
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: fortran
          Assignee: unassigned at gcc dot gnu.org
          Reporter: zeccav at gmail dot com

c gfortran 4.8.1 -ffpe-trap=zero or invalid produces SIGFPE
c Program received signal SIGFPE: Floating-point exception - erroneous
arithmetic operation.
c I have package glibc-2.17-4.fc19.x86_64
c under valgrind 3.8.1 it works
c must be compiled and run
      complex zero
      zero=(0e0,0e0)
      print *,zero**1e0
      end

[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0

[anlauf at gmx dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

Harald Anlauf <anlauf at gmx dot de> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |anlauf at gmx dot de

--- Comment #1 from Harald Anlauf <anlauf at gmx dot de> ---
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

[anlauf at gmx dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #3 from Harald Anlauf <anlauf at gmx dot de> ---
(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

[anlauf at gmx dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #5 from Harald Anlauf <anlauf at gmx dot de> ---
(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

[zeccav at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #7 from Vittorio Zecca <zeccav at gmail dot com> ---
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.00000000    , -0.00000000    ),
higher optimizations deliver (  0.00000000    ,  0.00000000    ).
Another mystery.

[Bug fortran/57749] -ffpe-trap=zero or invalid produces SIGFPE on complex zero ** 1e0

[dominiq at lps dot ens.fr]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #8 from Dominique d'Humieres <dominiq at lps dot ens.fr> ---
> 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.00000000    , -0.00000000    ),
> higher optimizations deliver (  0.00000000    ,  0.00000000    ).
> 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

[zeccav at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #9 from Vittorio Zecca <zeccav at gmail dot com> ---
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

[dominiq at lps dot ens.fr]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #10 from Dominique d'Humieres <dominiq at lps dot ens.fr> ---
> 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

[anlauf at gmx dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #11 from Harald Anlauf <anlauf at gmx dot de> ---
(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

[dominiq at lps dot ens.fr]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #13 from Dominique d'Humieres <dominiq at lps dot ens.fr> ---
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

[schwab@linux-m68k.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #14 from Andreas Schwab <schwab@linux-m68k.org> ---
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

[kargl at gcc dot gnu.org]

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

[zeccav at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #16 from Vittorio Zecca <zeccav at gmail dot com> ---
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

[sgk at troutmask dot apl.washington.edu]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #18 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Tue, Jul 02, 2013 at 07:21:54AM +0000, zeccav at gmail dot com wrote:
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749
> 
> --- Comment #16 from Vittorio Zecca <zeccav at gmail dot com> ---

> 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

[anlauf at gmx dot de]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #19 from Harald Anlauf <anlauf at gmx dot de> ---
(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

[zeccav at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #20 from Vittorio Zecca <zeccav at gmail dot com> ---
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

[zeccav at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57749

--- Comment #21 from Vittorio Zecca <zeccav at gmail dot com> ---
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.