-m64 changes the default from -mfp-math=387 to something else. That is clearly
a good idea, but has an unexpected downside for complex division, in that the
boundary case handling degrades badly. This was detected for C under 4.1.2 but
is almost certainly more general.
It is worth considering whether complex division should be a special case, and
use Intel arithmetic when it is available, irrespective of the arithmetic used
for most floating-point calculations. Sigh - another hack :-(
Take the following example:
#include <complex.h>
#include <stdio.h>
complex f(double complex a, double complex b) { return a/b;}
int main()
{ double complex X = 1.0e308;
int n;
for (n = 0; n<20; n++)
{ double complex a = X + X*I, b = X + (n/10.0)*X*I;
double complex r = f(a,b);
printf("(%.3e,%.3e)/(%.3e,%.3e) => (%.3e,%.3e)\n",a,b,r);
}
return 0;
}
This produces the following output with -m32 or -mfpmath=387:
(1.000e+308,1.000e+308)/(1.000e+308,0.000e+00) => (1.000e+00,1.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.000e+307) => (1.089e+00,8.911e-01)
(1.000e+308,1.000e+308)/(1.000e+308,2.000e+307) => (1.154e+00,7.692e-01)
(1.000e+308,1.000e+308)/(1.000e+308,3.000e+307) => (1.193e+00,6.422e-01)
(1.000e+308,1.000e+308)/(1.000e+308,4.000e+307) => (1.207e+00,5.172e-01)
(1.000e+308,1.000e+308)/(1.000e+308,5.000e+307) => (1.200e+00,4.000e-01)
(1.000e+308,1.000e+308)/(1.000e+308,6.000e+307) => (1.176e+00,2.941e-01)
(1.000e+308,1.000e+308)/(1.000e+308,7.000e+307) => (1.141e+00,2.013e-01)
(1.000e+308,1.000e+308)/(1.000e+308,8.000e+307) => (1.098e+00,1.220e-01)
(1.000e+308,1.000e+308)/(1.000e+308,9.000e+307) => (1.050e+00,5.525e-02)
(1.000e+308,1.000e+308)/(1.000e+308,1.000e+308) => (1.000e+00,0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.100e+308) => (9.502e-01,-4.525e-02)
(1.000e+308,1.000e+308)/(1.000e+308,1.200e+308) => (9.016e-01,-8.197e-02)
(1.000e+308,1.000e+308)/(1.000e+308,1.300e+308) => (8.550e-01,-1.115e-01)
(1.000e+308,1.000e+308)/(1.000e+308,1.400e+308) => (8.108e-01,-1.351e-01)
(1.000e+308,1.000e+308)/(1.000e+308,1.500e+308) => (7.692e-01,-1.538e-01)
(1.000e+308,1.000e+308)/(1.000e+308,1.600e+308) => (7.303e-01,-1.685e-01)
(1.000e+308,1.000e+308)/(1.000e+308,1.700e+308) => (6.941e-01,-1.799e-01)
(1.000e+308,1.000e+308)/(1.000e+308,inf) => (0.000e+00,-0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,inf) => (0.000e+00,-0.000e+00)
[ Actually, with -mfpmath=387, the last two lines are:
(1.000e+308,1.000e+308)/(nan,inf) => (nan,nan)
(1.000e+308,1.000e+308)/(nan,inf) => (nan,nan) ]
This produces the following output with -m64 or -mfpmath=sse:
(1.000e+308,1.000e+308)/(1.000e+308,0.000e+00) => (1.000e+00,1.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.000e+307) => (1.089e+00,8.911e-01)
(1.000e+308,1.000e+308)/(1.000e+308,2.000e+307) => (1.154e+00,7.692e-01)
(1.000e+308,1.000e+308)/(1.000e+308,3.000e+307) => (1.193e+00,6.422e-01)
(1.000e+308,1.000e+308)/(1.000e+308,4.000e+307) => (1.207e+00,5.172e-01)
(1.000e+308,1.000e+308)/(1.000e+308,5.000e+307) => (1.200e+00,4.000e-01)
(1.000e+308,1.000e+308)/(1.000e+308,6.000e+307) => (1.176e+00,2.941e-01)
(1.000e+308,1.000e+308)/(1.000e+308,7.000e+307) => (1.141e+00,2.013e-01)
(1.000e+308,1.000e+308)/(1.000e+308,8.000e+307) => (inf,1.220e-01)
(1.000e+308,1.000e+308)/(1.000e+308,9.000e+307) => (nan,0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.000e+308) => (nan,0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.100e+308) => (nan,-0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.200e+308) => (nan,-0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.300e+308) => (0.000e+00,-0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.400e+308) => (0.000e+00,-0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.500e+308) => (0.000e+00,-0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.600e+308) => (0.000e+00,-0.000e+00)
(1.000e+308,1.000e+308)/(1.000e+308,1.700e+308) => (0.000e+00,-0.000e+00)
(1.000e+308,1.000e+308)/(nan,inf) => (nan,nan)
(1.000e+308,1.000e+308)/(nan,inf) => (nan,nan)
Note the particularly horrible case, where a value that should be small
is actually infinite!
--
Summary: Poor behaviour of complex division in 64-bit on Intel
Product: gcc
Version: 4.1.2
Status: UNCONFIRMED
Severity: normal
Priority: P3
Component: c
AssignedTo: unassigned at gcc dot gnu dot org
ReportedBy: nmm1 at cam dot ac dot uk
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34071
