------- Comment #9 from sgk at troutmask dot apl dot washington dot edu
2009-06-01 14:56 -------
Subject: Re: Complex division by zero in gfortran returns wrong results
On Mon, Jun 01, 2009 at 08:35:05AM -0000, ghazi at gcc dot gnu dot org wrote:
>
>
> ------- Comment #2 from ghazi at gcc dot gnu dot org 2009-06-01 08:35 -------
> (In reply to comment #1)
> > Kaveh,
> > After looking into the problem, I think (nan + i nan) is
> > an acceptable result for z = (-0.1,-2.2)/(0.0,0.0)
> > because of the standard definition of complex division.
> > For z2 = (0.1,1)/0, I think that you are correct, and
> > gfortran should give (inf + i inf).
>
> Why is one different than the other? I don't know fortan promotion rules, but
> in C, the latter case promotes to (0.1,1.0)/(0.0,0.0) which looks very much
> like the first case. Does fortran handle type promotion differently?
>
> Regardless, I don't know that any "standard definition" of complex division
> applies here. The canonical algebraic formula is undefined mathematically in
> the presence of division by zero. So at least in C there are rules telling us
> what to do, and they say return Inf. Does fortran follow a standard here we
> can compare against or are we guessing? :-)
>
>
Let's deal with the z = (-0.1,-2.2)/(0.0,0.0) case. This program, IMHO, needs
to give a consistent answer.
subroutine sub(z1, z2)
implicit none
complex z1, z2
print *, z1 / z2
end subroutine sub
program a
implicit none
complex :: z1 = (-0.1,-2.2), z2 = (0.0,0.0)
complex :: z3 = (-0.1,-2.2) / (0.0,0.0)
call sub(z1, z2)
print *, z3
end program a
If one wants z3 to be inf or (inf + i inf) or anything other than (nan + i
nan),
then the complex division in sub() will cause an unacceptably large drop in
performance because GCC would need to generate code to check all the special
cases i.e.,
subroutine sub(z1, z2)
implicit none
complex z1, z2
if (real(z2) == 0. .and. aimag(z2) == 0.) then
if (real(z1) == 0. .and. aimag(z1) then
nan + i nan
else if (aimag(z1) == 0.) then
inf + i nan
else if (real(z1) == 0.) then
nan + i inf
end if
else
print *, z1 / z2
end if
end subroutine sub
As for the second case of z = (0.1, 1) / 0, Fortran indeed has promotion
rules (see Sec. 7.1.2 in Fortran 2003). This is transformed to
z = (0.1, 1) / (0., 0.), so once again my above argument for a consistent
result applies.
PS: Please do not consider -ffast-math as a method to disable all of the
checking. -ffast-math has too much baggage to be used with impunity.
If MPC returns inf or (inf + i inf) and the MPC developers do not consider
this to be a bug in their library, then gfortran will need to handle the
division by zero during constant folding as a special case.
--
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=40318
