[Cfp-interest 1472] complex multiplication and division

[paul zimmermann]

       Dear all,

I'm new to this list, I hope this is not off-topic.

I recently discovered that GCC (version 9.2.1) does not round correctly
complex multiplications or divisions:

$ cat complex.c
#include <stdio.h>
#include <complex.h>

int
main()
{
  double complex x, y, z;

  /* generated with ./mpcheck-double -seed 1473128 */
  x = -8.2287892462057963e+188 + I * 4.7082627796239037e+122;
  y = -1.0987942433164555e-323 + I * 6.4181344246145655e-324;
  printf ("x=(%.16e,%.16e)\n", creal (x), cimag (x));
  printf ("y=(%.16e,%.16e)\n", creal (y), cimag (y));
  z = x * y;
  printf ("x*y=(%.16e,%.16e)\n", creal (z), cimag (z));

  printf ("\n");

  x = -2.2994568145246667e-322 + I * -6.7148897446878843e-324;
  y = -2.9555380987538969e-207 + I * 2.1451122087941628e-255;
  printf ("x=(%.16e,%.16e)\n", creal (x), cimag (x));
  printf ("y=(%.16e,%.16e)\n", creal (y), cimag (y));
  z = x / y;
  printf ("x/y=(%.16e,%.16e)\n", creal (z), cimag (z));
}

$ gcc complex.c 
$ ./a.out 
x=(-8.2287892462057963e+188,4.7082627796239037e+122)
y=(-9.8813129168249309e-324,4.9406564584124654e-324)
x*y=(8.1311241468363422e-135,-4.0655620734181711e-135)

x=(-2.3221085354538588e-322,-4.9406564584124654e-324)
y=(-2.9555380987538969e-207,2.1451122087941628e-255)
x/y=(7.8568046083821340e-116,1.6716605549749219e-117)

The correct results (as computed with MPFR) should be:

x*y=9.04174625319528e-135 - 5.28134755339716e-135*I
x/y=7.78016299466468e-116 + 2.27196859601268e-117*I

i.e., even the first digits are wrong!

My question is thus: is there any correct rounding requirement for
complex operations, as there is for floating-point ones?

Best regards,
Paul Zimmermann

PS: with ICC 19.0.4.243 (gcc version 9.2.0 compatibility) I get:

$ icc complex.c; ./a.out 
x=(-8.2287892462057963e+188,4.7082627796239037e+122)
y=(-9.8813129168249309e-324,4.9406564584124654e-324)
x*y=(0.0000000000000000e+00,-0.0000000000000000e+00)

x=(-2.3221085354538588e-322,-4.9406564584124654e-324)
y=(-2.9555380987538969e-207,2.1451122087941628e-255)
x/y=(7.8568046083821340e-116,1.6716605549749219e-117)