This is the mail archive of the gcc-bugs@gcc.gnu.org mailing list for the GCC project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]