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

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

[ghazi at gcc dot gnu dot org]

Complex division by zero in gfortran returns NaN.  This is expected for 0/0,
but other finite/zero should return Inf.  This impacts the testcase
gfortran.dg/real_const_3.f90 in two values incorrectly computed:

  complex :: z = (-0.1,-2.2)/(0.0,0.0)
  complex :: z2 = (0.1,1)/0

See: http://gcc.gnu.org/ml/fortran/2009-05/msg00423.html

This should be fixed in gcc-4.5 by using MPC for division, but older versions
of GCC should add special case handling in the fortran frontend simplification
code.


-- 
           Summary: Complex division by zero in gfortran returns wrong
                    results
           Product: gcc
           Version: 4.3.4
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: fortran
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: ghazi at gcc dot gnu dot org


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

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

[kargl at gcc dot gnu dot org]

------- Comment #1 from kargl at gcc dot gnu dot org  2009-06-01 06:54 -------
(In reply to comment #0)
> Complex division by zero in gfortran returns NaN.  This is expected for 0/0,
> but other finite/zero should return Inf.  This impacts the testcase
> gfortran.dg/real_const_3.f90 in two values incorrectly computed:
> 
>   complex :: z = (-0.1,-2.2)/(0.0,0.0)
>   complex :: z2 = (0.1,1)/0
> 
> See: http://gcc.gnu.org/ml/fortran/2009-05/msg00423.html
> 
> This should be fixed in gcc-4.5 by using MPC for division, but older versions
> of GCC should add special case handling in the fortran frontend simplification
> code.
> 

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


-- 


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

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

[ghazi at gcc dot gnu dot org]

------- 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? :-)


-- 


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

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

[dominiq at lps dot ens dot fr]

------- Comment #3 from dominiq at lps dot ens dot fr  2009-06-01 10:19 -------
In http://gcc.gnu.org/ml/fortran/2009-06/msg00006.html, Dennis Wassel 
wrote:
> Complex numbers have a well-defined concept of infinity, which I like to
> visualise as the "infinite-diameter" ring around the finite complex
> numbers, so its argument is, of course, not well-defined.  This is
> conceptually somewhat different from the two points +/-Inf on the real
> line, and it is not signed.

My understanding of infinity in the complex plane is what is called (I
call?)  "directed inifinity": if abs((a,b)) goes to +Inf and atan2(a,b) has
a defined value in this limit, then (a,b) goes to infinity in the direction
given by atan2(a,b).  However (+/-Inf,+/-Inf) defines only four directions
and is unable to represent general "directed inifinity".  So I think that
from a "mathematical" point of view the problem is ill-posed,
(+/-Inf,+/-Inf) is in the same class of "undefined" values as 0/0 or
Inf/Inf, and should give NaNs.

Now atan2(a,b) is DEFINED such that atan2(+0,+/-0)=+/-0,
atan2(-0,+/-0)=+/-Pi, atan2(+Inf,+/-Inf)=+/-Pi/4,
atan2(-Inf,+/-Inf)=+/-3*Pi/4 (it seems that it is even built in the Intel
hardware).  With this definition of atan2, it is possible to give a
definition of (+/-Inf,+/-Inf) as the directions of the corners of the
"infinite square".  I have nothing against this defintion, except it should
be documented.

Part of the problem originates from the optimization of (a,b)/(c,0.0) as
(a/c,b/c), see pr30482, and it limit when c goes to 0 (the latter giving
(+/-Inf,+/-Inf) if b/=0.0 or has a unknown value at compile time,
(+/-Inf,+/-0) if b known to be zero at compile time, optimized as
(a/c,+/-0.0), or (+/-Inf,-/+Nan) if b==0 at run time).

Note that due to the definition of atan2(+/-0,+/-0),
(+/-0,+/-0)=0*exp((+/-k*Pi,+/-0)) (k being 0 or 1), then to be consistent
1.0/(+/-0.0,+/-0.0) should give 1/0.0*exp((-/+k*Pi,-/+0))=(+/-Inf,-/+Nan).

Last remark, if I remember correctly (I cannot find the right pointer) C99 
defines precisely how complex multiplications should behave in the 
exception limits, with the drawback of a large increase of the computation 
cost avoided with the -fcx-fortran-rules option.

As far as I can tell:

(1) Without the IEEE module, using or producing Inf or NaN makes the code 
non (fortran) standard comforming, so the "processor" can give any answer.

(2) With the IEEE module, Inf and NaN are part of the numerical model, 
however I did not find any definition of the values that intrinsincs 
should return in such cases.

To conclude, I'll repeat that in my opinion the only relevant answer for
overflow, divide by zero, or invalid argument exceptions is "call abort()".
Since I know that it is a lost battle, I think at least the behavior shall
not increase the computation time with complex numbers (or provide a way to
avoid the penalty).  Now if the fxxxx standard does not define the expected
result, but the C99 does, the best answer to this ill-posed problem is
probably to follow the C99 standard if possible (I have no idea of what
does mpc).  In this case the results at compile and run times should be the
same ("least surprising approach").


-- 


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

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

[burnus at gcc dot gnu dot org]

------- Comment #4 from burnus at gcc dot gnu dot org  2009-06-01 11:09 -------
Regarding the run-time evaluation, note that Fortran sets (internally) the
-fcx-fortran-rules flag, which modifies the complex evaluation.


-- 

burnus at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |burnus at gcc dot gnu dot
                   |                            |org
           Keywords|                            |wrong-code


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

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

[dominiq at lps dot ens dot fr]

------- Comment #5 from dominiq at lps dot ens dot fr  2009-06-01 11:24 -------
(In reply to comment #2)
> Does fortran follow a standard here we can compare against 
> or are we guessing? :-)

What the fortran standard says is "you shall not divide by zero"! anything else
is just a matter of choice of the implementation, but it won't make the code
valid.


-- 


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

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

[dennis dot wassel at googlemail dot com]

------- Comment #6 from dennis dot wassel at googlemail dot com  2009-06-01 12:14 -------
(In reply to comment #3)

> My understanding of infinity in the complex plane is what is called (I
> call?)  "directed inifinity": if abs((a,b)) goes to +Inf and atan2(a,b) has
> a defined value in this limit, then (a,b) goes to infinity in the direction
> given by atan2(a,b).  However (+/-Inf,+/-Inf) defines only four directions
> and is unable to represent general "directed inifinity".  So I think that
> from a "mathematical" point of view the problem is ill-posed,
> (+/-Inf,+/-Inf) is in the same class of "undefined" values as 0/0 or
> Inf/Inf, and should give NaNs.

Then this is the gist of the matter - my FA textbook does not require the
argument to converge, but just the modulus, so our understandings of infinity
differ.
I agree with Dominique in that mathematically (+/-Inf,+/-Inf) is not
well-defined, because the limit of a/b is unknown.

> Now atan2(a,b) is DEFINED such that atan2(+0,+/-0)=+/-0,
> atan2(-0,+/-0)=+/-Pi, atan2(+Inf,+/-Inf)=+/-Pi/4,
> atan2(-Inf,+/-Inf)=+/-3*Pi/4 (it seems that it is even built in the Intel
> hardware).  With this definition of atan2, it is possible to give a
> definition of (+/-Inf,+/-Inf) as the directions of the corners of the
> "infinite square".  I have nothing against this defintion, except it should
> be documented.

The mathematicians among us will writhe in agony, but I think this serves the
purpose best :)

> [...]
> In this case the results at compile and run times should be the
> same ("least surprising approach").

I strongly support that point, because everything else would be a very nasty
gotcha for an average Fortran user like myself.


-- 


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

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

[dominiq at lps dot ens dot fr]

------- Comment #7 from dominiq at lps dot ens dot fr  2009-06-01 12:25 -------
> Then this is the gist of the matter - my FA textbook does not require the
> argument to converge, but just the modulus, so our understandings of infinity differ.

Think of something like \rho\exp(i\rho\sin(\pi\rho)): it is going to infinity
when \rho goes to +infinity, but you cannot define a directed infinity since
the argument keeps oscillating when \rho increases.


-- 


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

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

[jb at gcc dot gnu dot org]

------- Comment #8 from jb at gcc dot gnu dot org  2009-06-01 14:21 -------
Whatever you do, as long as the Fortran standard is silent on this matter, can
you hide it behind some -fC99-wankery or such option and not enable it by
default, so that those of us who care less about which of (NaN, NaN) and (Inf,
Inf) is the better error indicator than the roughly factor of 4 in performance
for complex arithmetic can get on with our lives?


-- 

jb at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jb at gcc dot gnu dot org


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

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

[sgk at troutmask dot apl dot washington dot edu]

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

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

[ghazi at gcc dot gnu dot org]

------- Comment #10 from ghazi at gcc dot gnu dot org  2009-06-01 18:14 -------
(In reply to comment #9)
> 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.

I believe the goals for MPC are to follow C99 rules for special cases.  Thus
the return value of (inf + i inf) is intentional for MPC and not a bug in their
thinking.

I entirely agree that the compile-time and runtime results should be identical.
 If it is your intention to preserve the existing runtime behavior, then we
should do the same in the fortran folder and special case this if/when
converting complex division to use MPC.

Does this mean this PR should be closed as "invalid" ?


-- 


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

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

[sgk at troutmask dot apl dot washington dot edu]

------- Comment #11 from sgk at troutmask dot apl dot washington dot edu  2009-06-01 19:16 -------
Subject: Re:  Complex division by zero in gfortran returns wrong results

On Mon, Jun 01, 2009 at 06:14:52PM -0000, ghazi at gcc dot gnu dot org wrote:
> 
> ----- Comment #10 from ghazi at gcc dot gnu dot org  2009-06-01 18:14 -------
> (In reply to comment #9)
> > 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.
> 
> I believe the goals for MPC are to follow C99 rules for special cases.
> Thus the return value of (inf + i inf) is intentional for MPC and not
> a bug in their thinking.
> 
> I entirely agree that the compile-time and runtime results should be
> identical.  If it is your intention to preserve the existing runtime
> behavior, then we should do the same in the fortran folder and special
> case this if/when converting complex division to use MPC.

Oh yuck.  I just checked n1124.pdf  In Annex G.5.1, it explicitly
defines this behavior:
   "if the first operand is a nonzero finite number or an infinity
    and the second operand is a zero, then the result of the / operator
    is an infinity."

Combining this with G.3:
   "A complex or imaginary value is a zero if each of its parts is a zero."

What is disturbing is Example 2 in G.5.1 on page 470!  Does gcc's runtime
implementation of complex division mirror Example 2?  I can understand
the need to avoid under/overflow, but _Cdivd() seems overly complicated.   

> Does this mean this PR should be closed as "invalid" ?

I think we should leave it open, perhaps in a suspended state, 
as a reminder (to me) that soemthing needs to be done with
gfortran an z/(0,0).


-- 


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

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

[burnus at gcc dot gnu dot org]

------- Comment #12 from burnus at gcc dot gnu dot org  2009-06-01 20:33 -------
> Oh yuck.  I just checked n1124.pdf  In Annex G.5.1, it explicitly
> defines this behavior:

Note: Annex G is "only" informative and not normative; still it makes probably
sense to follow the informative parts of a standard unless they clash with a
normative part - or are intolerably expensive in implementation, compile, or
run time.


-- 


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

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

[ghazi at gcc dot gnu dot org]

------- Comment #13 from ghazi at gcc dot gnu dot org  2009-06-02 15:16 -------
(In reply to comment #11)
> What is disturbing is Example 2 in G.5.1 on page 470!  Does gcc's runtime
> implementation of complex division mirror Example 2?  I can understand
> the need to avoid under/overflow, but _Cdivd() seems overly complicated.   

Here is GCC's runtime implementation of complex division from libgcc2.c.  It
looks like it does mirror example 2.  While the runtime evaluation seems to be
fine, the middle-end folder still has bugs.  See PR30789.

------------------------------------------------------------

#if defined(L_divsc3) || defined(L_divdc3) \
    || defined(L_divxc3) || defined(L_divtc3)

CTYPE
CONCAT3(__div,MODE,3) (MTYPE a, MTYPE b, MTYPE c, MTYPE d)
{
  MTYPE denom, ratio, x, y;
  CTYPE res;

  /* ??? We can get better behavior from logarithmic scaling instead of
     the division.  But that would mean starting to link libgcc against
     libm.  We could implement something akin to ldexp/frexp as gcc builtins
     fairly easily...  */
  if (FABS (c) < FABS (d))
    {
      ratio = c / d;
      denom = (c * ratio) + d;
      x = ((a * ratio) + b) / denom;
      y = ((b * ratio) - a) / denom;
    }
  else
    {
      ratio = d / c;
      denom = (d * ratio) + c;
      x = ((b * ratio) + a) / denom;
      y = (b - (a * ratio)) / denom;
    }

  /* Recover infinities and zeros that computed as NaN+iNaN; the only cases
     are nonzero/zero, infinite/finite, and finite/infinite.  */
  if (isnan (x) && isnan (y))
    {
      if (c == 0.0 && d == 0.0 && (!isnan (a) || !isnan (b)))
        {
          x = COPYSIGN (INFINITY, c) * a;
          y = COPYSIGN (INFINITY, c) * b;
        }
      else if ((isinf (a) || isinf (b)) && isfinite (c) && isfinite (d))
        {
          a = COPYSIGN (isinf (a) ? 1 : 0, a);
          b = COPYSIGN (isinf (b) ? 1 : 0, b);
          x = INFINITY * (a * c + b * d);
          y = INFINITY * (b * c - a * d);
        }
      else if ((isinf (c) || isinf (d)) && isfinite (a) && isfinite (b))
        {
          c = COPYSIGN (isinf (c) ? 1 : 0, c);
          d = COPYSIGN (isinf (d) ? 1 : 0, d);
          x = 0.0 * (a * c + b * d);
          y = 0.0 * (b * c - a * d);
        }
    }

  __real__ res = x;
  __imag__ res = y;
  return res;
}
#endif /* complex divide */


-- 


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

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

[pault at gcc dot gnu dot org]

------- Comment #14 from pault at gcc dot gnu dot org  2009-12-18 14:42 -------
Hi guys!

Do you want to suspend this PR or to junk it?

Let's get it out of the unconfirmed list.

Thanks

Paul


-- 


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

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

[sgk at troutmask dot apl dot washington dot edu]

------- Comment #15 from sgk at troutmask dot apl dot washington dot edu  2009-12-18 15:52 -------
Subject: Re:  Complex division by zero in gfortran returns wrong results

On Fri, Dec 18, 2009 at 02:42:15PM -0000, pault at gcc dot gnu dot org wrote:
> 
> Do you want to suspend this PR or to junk it?
> 
> Let's get it out of the unconfirmed list.
> 

Now that MPC is required by gcc, I'll take a look
at making gfortran give a consistent result when
comparing its constant folding with generated code.


-- 


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

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

[ghazi at gcc dot gnu dot org]

------- Comment #16 from ghazi at gcc dot gnu dot org  2009-12-18 17:16 -------
(In reply to comment #15)
> Now that MPC is required by gcc, I'll take a look
> at making gfortran give a consistent result when
> comparing its constant folding with generated code.

BTW, I put in some special-case code in arith.c:gfc_arith_divide() to ensure we
still got the old fortran divide behavior.  You probably want to remove it if
you're working on this PR.  See note #3 from this patch:

http://gcc.gnu.org/ml/gcc-patches/2009-06/msg01368.html


-- 


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

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

[kargl at gcc dot gnu dot org]

------- Comment #17 from kargl at gcc dot gnu dot org  2009-12-19 19:45 -------
Closing as won't fix.  Currently, gfortran produces nan + i nan for
finite-complex / 0 or finite-complex / (0,0) with its FE constanting
folding and when the middle-end generates code.  The Fortran standard
does not provide any guidance with respect to exceptional FP values,
and trying to get the same output of an equivalent C program is not
required by either standard.  


-- 

kargl at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |RESOLVED
         Resolution|                            |WONTFIX


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