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