[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