[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.