complex division is inaccurate

Bug #1266571 reported by Michael Hudson-Doyle
6
This bug affects 1 person
Affects Status Importance Assigned to Milestone
Linaro GCC
Won't Fix
Undecided
Unassigned
gcc
In Progress
Medium

Bug Description

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).

Revision history for this message
Michael Hudson-Doyle (mwhudson) wrote :

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).

Revision history for this message
Matthew Gretton-Dann (matthew-gretton-dann) wrote :

Michael,

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

Thanks,

Matt

Revision history for this message
In , Michael-hudson (michael-hudson) wrote :

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).

Revision history for this message
Michael Hudson-Doyle (mwhudson) wrote :
Revision history for this message
In , Pinskia (pinskia) wrote :
Revision history for this message
In , Michael-hudson (michael-hudson) wrote :

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...

Revision history for this message
In , Joseph-codesourcery (joseph-codesourcery) wrote :

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).

Revision history for this message
In , Michael-hudson (michael-hudson) wrote :

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.

Changed in gcc-linaro:
status: New → Won't Fix
Revision history for this message
In , Pinskia (pinskia) wrote :

Also see PR 19974.

Changed in gcc:
importance: Unknown → Medium
status: Unknown → New
Revision history for this message
In , Wilco-z (wilco-z) wrote :

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.

Revision history for this message
In , Wilco-z (wilco-z) wrote :

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)

Revision history for this message
In , Joseph-codesourcery (joseph-codesourcery) wrote :

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).

Changed in gcc:
status: New → Confirmed
Changed in gcc:
status: Confirmed → In Progress
Revision history for this message
In , Patrick-mcgehearty (patrick-mcgehearty) wrote :

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.

Revision history for this message
In , Pinskia (pinskia) wrote :

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

Just making it easier for links here.

Revision history for this message
In , Patrick-mcgehearty (patrick-mcgehearty) wrote :

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.
>

To post a comment you must log in.
This report contains Public information  
Everyone can see this information.

Other bug subscribers

Remote bug watches

Bug watches keep track of this bug in other bug trackers.