The constant folding of arithmetic on complex numbers is inexact and doesn't
support any of the C99 special cases and avoiding overflow done in libgcc.
Compile and run the following testcase with -std=c99.
#include <float.h>
int printf (const char *, ...);
#define C1 (DBL_MAX * 0.5 + DBL_MAX * 0.5i)
#define C2 (DBL_MAX * 0.25 + DBL_MAX * 0.25i)
_Complex double c1 = C1;
_Complex double c2 = C2;
_Complex double c = C1 / C2;
int
main (void)
{
_Complex double cr = c1 / c2;
printf ("%f %f\n%f %f\n", __real__ c, __imag__ c, __real__ cr, __imag__ cr);
return 0;
}
It prints:
nan nan
2.000000 0.000000
That is, the libgcc code handled the division correctly, the folding allowed it
to overflow.
It should be reasonably straightforward to fold complex multiplication exactly
using MPFR: compute ac and bd with results in double the input precision, then
compute ac-bd with result in the original input precision (so rounding just
once to that precision) and similarly with the imaginary part. Getting
division exact (correctly rounded) may be trickier, but doing better than at
present should be straightforward (simply by taking advantage of the greater
exponent range in MPFR). The folding code also needs to handle the C99 Annex G
special cases like the libgcc code does.
--
Summary: complex folding inexact
Product: gcc
Version: 4.3.0
Status: UNCONFIRMED
Severity: normal
Priority: P3
Component: middle-end
AssignedTo: unassigned at gcc dot gnu dot org
ReportedBy: jsm28 at gcc dot gnu dot org
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=30789
