[Bug c/34071] Poor behaviour of complex division in 64-bit on Intel

[rguenth at gcc dot gnu dot org]

------- Comment #1 from rguenth at gcc dot gnu dot org  2007-11-12 13:41 -------
We open-code complex division in libgcc and benefit from the extra precision on
x87 in this case (so this is related to PR323).  But surely __divdc3 has
QOI issues here:

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

  /* ??? 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);
        }
    }

  return x + I * y;
}
#endif /* complex divide */


I bet we have a duplicate bug of this somewhere, but I can't find it right
now.


-- 


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