[Bug fortran/40318] Complex division by zero in gfortran returns wrong results

["sgk at troutmask dot apl dot washington dot edu" <gcc-bugzilla at gcc dot gnu dot org>]

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