This is the mail archive of the
gcc-bugs@gcc.gnu.org
mailing list for the GCC project.
[Bug fortran/40318] Complex division by zero in gfortran returns wrong results
- From: "ghazi at gcc dot gnu dot org" <gcc-bugzilla at gcc dot gnu dot org>
- To: gcc-bugs at gcc dot gnu dot org
- Date: 2 Jun 2009 15:16:38 -0000
- Subject: [Bug fortran/40318] Complex division by zero in gfortran returns wrong results
- References: <bug-40318-578@http.gcc.gnu.org/bugzilla/>
- Reply-to: gcc-bugzilla 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