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

["nmm1 at cam dot ac dot uk" <gcc-bugzilla at gcc dot gnu dot org>]

-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