complex division is inaccurate

[Michael Hudson-Doyle]

Hi, it seems details in the rounding of complex division are wrong:

ubuntu@arm64:~$ cat cplx.c
#include <stdio.h>

int main(int argc, char** argv)
{
 __complex double c = 1.0 + 3.0i;
 printf("%g\n", __imag__ (c/c));
}
ubuntu@arm64:~$ gcc cplx.c -o cplx
ubuntu@arm64:~$ ./cplx
-1.66533e-17

This happens with 4.8 and 4.9 (as of a couple of weeks ago).

[Michael Hudson-Doyle (mwhudson)]

This is because libgcc2.c is compiled in a way that lets the compiler used fused multiply add instructions. It shouldn't be!

(in this particular case, it's because if x is the closest float64 to 1/3, evaluating "3*x - 1" using fnmsub does not yield 0 because of the lack of intermediate rounding. Evaluating 3*x, rounding it, and then subtracting 1 does yield 0).

[Matthew Gretton-Dann (matthew-gretton-dann)]

Michael,

This is a generic issue can you please raise this in GCC BugZilla: http://gcc.gnu.org/bugzilla?

Thanks,

Matt

[Michael-hudson (michael-hudson)]

Hi, it seems details in the rounding of complex division are wrong:

ubuntu@arm64:~$ cat cplx.c
#include <stdio.h>

int main(int argc, char** argv)
{
 __complex double c = 1.0 + 3.0i;
 printf("%g\n", __imag__ (c/c));
}
ubuntu@arm64:~$ gcc cplx.c -o cplx
ubuntu@arm64:~$ ./cplx
-1.66533e-17

This is because libgcc2.c is compiled in a way that lets the compiler used fused multiply add instructions. It shouldn't be!

(in this particular case, it's because if x is the closest float64 to 1/3, evaluating "3*x - 1" using fnmsub does not yield 0 because of the lack of intermediate rounding. Evaluating 3*x, rounding it, and then subtracting 1 does yield 0).

[Michael Hudson-Doyle (mwhudson)]

Sure: http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714

[Pinskia (pinskia)]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=19974#c6

[Michael-hudson (michael-hudson)]

Oh right. I saw that bug but didn't read enough of the comments to see that the same thing was mentioned.

I still think fma is inappropriate for this function...

[Joseph-codesourcery (joseph-codesourcery)]

You'll find other cases where fused operations make complex multiplication
and division *more* accurate.

There is, I suppose, a case for having a version of complex multiplication
and division that follow the same accuracy goals as glibc libm (i.e. real
and imaginary parts within a few ulp of the mathematically exact values,
with corresponding exceptions, avoiding spurious overflow / underflow) -
or even correctly rounding versions. But (a) people wouldn't want them as
default, on grounds of speed, and (b) they would be tricky to implement in
libgcc because of needing to change the rounding mode temporarily (libgcc
can't depend on libm) and because of depending on details of the
floating-point formats that libgcc otherwise doesn't need to depend on.
Maybe people interested in complex arithmetic should propose new
facilities for accurate multiplication / division when the C standard is
next being revised (there's some interest in other new functions, see
N1704, but no current project in which they could be added to the
standard).

[Michael-hudson (michael-hudson)]

OK, so I should stop talking about accurate and instead talk about surprising :-)

I think it is pretty surprising that x/x != 1+0i for (I think) all x where the ratio between the real and imaginary parts is not a power of 2 (put another way, an error of 1ulp when the accurate result is 0 is more surprising than when the accurate result is non-zero).

In general, fma makes it hard to reason about expressions like "a*b-c*d" -- which of the multiplications is being done at the higher precision? But I guess arguing to fp-contract to default to off is a war I don't want to get into fighting.

[Pinskia (pinskia)]

Also see PR 19974.

[Wilco-z (wilco-z)]

The question is whether the algorithm used in __divdc3 is accurate - it appears to want to use FMA explictly, otherwise you'd see uses of TRUNC like in __muldc3. When I build the example with -Ofast (adding volatile etc to avoid it optimizing away), I get:

 fmul d3, d3, d2
 fmul d2, d2, d2
 fnmsub d0, d0, d1, d3
 fmadd d1, d1, d1, d2
 fdiv d0, d0, d1

> ./a.out
0

With -O3:

 bl __divdc3

> ./a.out
-1.66533e-17

So this bit of the code needs to be investigated for accuracy when using FMA:

  if (FABS (c) < FABS (d))
    {
      ratio = c / d;
      denom = (c * ratio) + d;
      x = ((a * ratio) + b) / denom;
      y = ((b * ratio) - a) / denom;
    }
  else
    {
      ratio = d / c;
      denom = (d * ratio) + c;
      x = ((b * ratio) + a) / denom;
      y = (b - (a * ratio)) / denom;
    }

As Richard already pointed out, this is a generic issue.

[Wilco-z (wilco-z)]

Btw this is also totally broken in libgcc2.c:

#define isnan(x) __builtin_expect ((x) != (x), 0)
#define isfinite(x) __builtin_expect (!isnan((x) - (x)), 1)
#define isinf(x) __builtin_expect (!isnan(x) & !isfinite(x), 0)

[Joseph-codesourcery (joseph-codesourcery)]

The algorithm isn't expecting to use FMA; it's just treating
floating-point numbers as approximate real numbers. I'm not sure why
multiplication has TRUNC, but my guess is it's about excess range not
precision - excess precision would improve results, but excess range in
only one of two numbers added/subtracted could result in an infinity of
the wrong sign as a result on overflow.

The division for scaling is, as commented, inferior to scalbn / logb
approaches as in Annex G, but to use such logarithmic scaling would
require appropriate built-in functions for scalbn / logb supported with
inline expansion on the target. Even then, accuracy in the sense of at
most a few ulps error in real and imaginary parts separately requires fmma
(a*b + c*d) as the atomic operation, not fma; without that, you have the
more complicated task (if you want that accuracy) of constructing
(approximate) fmma out of available operations such as fma.

It's true that now we have full __builtin_isnan / __builtin_isfinite /
__builtin_isinf those macros for the complex functions should probably use
them.

As noted above, one could consider proposing complex multiplication /
division with better / better defined accuracy (or other new facilities)
for the next revision of the C standard. Proposals for C2x can be made
over the next few years (we don't actually have a C2x working draft yet,
just a C17 bug-fix version being reviewed, but you can still make C2x
proposals).

[Patrick-mcgehearty (patrick-mcgehearty)]

The patch 54f0224d55a1b56dde092460ddf76913670e6efc
"Practical improvement to libgcc complex divide"
made substantial improvements to the algorithms
used and thus the accuracy of complex divide for
IEEE-754 precisions for half, float, double, and
long-double precision.

In particular, when examining randomly generated
values uniformly distributed over the full exponent
range, 2 bit errors dropped from 1.88% to less than
one in 10,000. The particular input values discussed
in this problem report no longer show any errors
whether or not FMA is used.

See the patch for more complete details.

[Pinskia (pinskia)]

(In reply to Patrick McGehearty from comment #9)
r12-228-g54f0224d55a1b56dde092460ddf76913670e6efc

Just making it easier for links here.

[Patrick-mcgehearty (patrick-mcgehearty)]

Thanks!

On 9/9/2021 5:54 PM, pinskia at gcc dot gnu.org wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714
>
> --- Comment #10 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
> (In reply to Patrick McGehearty from comment #9)
> r12-228-g54f0224d55a1b56dde092460ddf76913670e6efc
>
> Just making it easier for links here.
>