[Cfp-interest 1474] Re: complex multiplication and division

[Jim Thomas]

Paul,

C has no requirement or recommendation for correctly rounding complex multiply or divide, even with Annex G support. Annex G support includes correct rounding for complex addition and subtraction.

- Jim Thomas

> On Feb 7, 2020, at 1:39 AM, paul zimmermann <Paul.Zimmermann at inria.fr> wrote:
> 
>       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)
> _______________________________________________
> Cfp-interest mailing list
> Cfp-interest at oakapple.net
> http://mailman.oakapple.net/mailman/listinfo/cfp-interest