[Bug middle-end/30789] New: complex folding inexact

[Bug middle-end/30789] New: complex folding inexact

[jsm28 at gcc dot gnu dot org]

The constant folding of arithmetic on complex numbers is inexact and doesn't
support any of the C99 special cases and avoiding overflow done in libgcc.

Compile and run the following testcase with -std=c99.

#include <float.h>

int printf (const char *, ...);

#define C1 (DBL_MAX * 0.5 + DBL_MAX * 0.5i)
#define C2 (DBL_MAX * 0.25 + DBL_MAX * 0.25i)

_Complex double c1 = C1;
_Complex double c2 = C2;
_Complex double c = C1 / C2;

int
main (void)
{
  _Complex double cr = c1 / c2;
  printf ("%f %f\n%f %f\n", __real__ c, __imag__ c, __real__ cr, __imag__ cr);
  return 0;
}

It prints:

nan nan
2.000000 0.000000

That is, the libgcc code handled the division correctly, the folding allowed it
to overflow.

It should be reasonably straightforward to fold complex multiplication exactly
using MPFR: compute ac and bd with results in double the input precision, then
compute ac-bd with result in the original input precision (so rounding just
once to that precision) and similarly with the imaginary part.  Getting
division exact (correctly rounded) may be trickier, but doing better than at
present should be straightforward (simply by taking advantage of the greater
exponent range in MPFR).  The folding code also needs to handle the C99 Annex G
special cases like the libgcc code does.


-- 
           Summary: complex folding inexact
           Product: gcc
           Version: 4.3.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: middle-end
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: jsm28 at gcc dot gnu dot org


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

[Bug middle-end/30789] complex folding inexact

[steven at gcc dot gnu dot org]

-- 

steven at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
OtherBugsDependingO|                            |27676
              nThis|                            |
             Status|UNCONFIRMED                 |NEW
     Ever Confirmed|0                           |1
   Last reconfirmed|0000-00-00 00:00:00         |2008-06-01 21:55:09
               date|                            |


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

[Bug middle-end/30789] complex folding inexact

[ghazi at gcc dot gnu dot org]

------- Comment #1 from ghazi at gcc dot gnu dot org  2009-07-21 17:20 -------
Joseph - I'm working on this one, but I'd appreciate it if you could help
compile a list of good test inputs beyond the one in the first comment.  I.e.
especially for the annex G stuff.  That way I can be more confident I'm writing
it correctly.

You can just list the inputs, you don't need to create an actual testcase. 
I'll do that.  Thanks!


-- 

ghazi at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
         AssignedTo|unassigned at gcc dot gnu   |ghazi at gcc dot gnu dot org
                   |dot org                     |
             Status|NEW                         |ASSIGNED
   Last reconfirmed|2008-06-01 21:55:09         |2009-07-21 17:20:48
               date|                            |


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

[Bug middle-end/30789] complex folding inexact

[joseph at codesourcery dot com]

------- Comment #2 from joseph at codesourcery dot com  2009-07-26 17:32 -------
Subject: Re:  complex folding inexact

The example in this bug deals with excess overflow for division.  For 
infinities computing as NaN + iNaN, an example is (NaN + iInf) * (NaN 
+iInf) (where NaN +iInf is obtained as (__builtin_inf () * (0.0 + 1.0i))) 
- this works if computed at runtime but not constant folded.  ("Works" 
here means that at least one part of the result is an infinity - I do not 
believe there is a specific requirement beyond that.)

For division producing NaN + iNaN, (1.0 + 0.0i) / (0.0 + 0.0i) should 
produce an infinity, (__builtin_inf () * (0.0 + 1.0i)) / (1.0 + 0.0i) 
should produce an infinity, (1.0 + 0.0i) / (__builtin_inf () * (0.0 + 
1.0i)) should produce a zero (in which each part may have either sign).

All the above can usefully be tested for both constant folding and runtime 
evaluation.

There are also cases for exact rounding where you'd expect MPC to produce 
the right results but would *not* expect operations executed at runtime to 
produce exactly those results.  For example, (1.0 + DBL_EPSILON + 1.0i) * 
(1.0 - DBL_EPSILON + 1.0i), which would only work at runtime if the code 
happens to use exactly the right fused multiply-add operation.


-- 


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

[Bug middle-end/30789] complex folding inexact

[ghazi at gcc dot gnu dot org]

------- Comment #3 from ghazi at gcc dot gnu dot org  2009-08-12 22:28 -------
(In reply to comment #2)

Joseph - Thanks for your reply and testvalues.

> There are also cases for exact rounding where you'd expect MPC to produce 
> the right results but would *not* expect operations executed at runtime to 
> produce exactly those results.  For example, (1.0 + DBL_EPSILON + 1.0i) * 
> (1.0 - DBL_EPSILON + 1.0i), which would only work at runtime if the code 
> happens to use exactly the right fused multiply-add operation.

What is the "right result" for this case?  GCC with MPC produces:
-4.93038065763132378382330353301741393545754021943139377981e-32 + 2.0i)

Unpatched GCC as well as the runtime on my x86_64 box says:
0.0 + 2.0i

So the runtime here is not using the fused instruction?

Is the MPC value correct?

Thanks.


-- 


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

[Bug middle-end/30789] complex folding inexact

[joseph at codesourcery dot com]

------- Comment #4 from joseph at codesourcery dot com  2009-08-13 01:25 -------
Subject: Re:  complex folding inexact

On Wed, 12 Aug 2009, ghazi at gcc dot gnu dot org wrote:

> 
> 
> ------- Comment #3 from ghazi at gcc dot gnu dot org  2009-08-12 22:28 -------
> (In reply to comment #2)
> 
> Joseph - Thanks for your reply and testvalues.
> 
> > There are also cases for exact rounding where you'd expect MPC to produce 
> > the right results but would *not* expect operations executed at runtime to 
> > produce exactly those results.  For example, (1.0 + DBL_EPSILON + 1.0i) * 
> > (1.0 - DBL_EPSILON + 1.0i), which would only work at runtime if the code 
> > happens to use exactly the right fused multiply-add operation.
> 
> What is the "right result" for this case?  GCC with MPC produces:
> -4.93038065763132378382330353301741393545754021943139377981e-32 + 2.0i)
> 
> Unpatched GCC as well as the runtime on my x86_64 box says:
> 0.0 + 2.0i
> 
> So the runtime here is not using the fused instruction?
> 
> Is the MPC value correct?

It looks correct.  You expect real part -DBL_EPSILON*DBL_EPSILON, 
imaginary part 2.0.

I believe existing x86_64 hardware does not have fused multiply-add 
instructions.


-- 


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

[Bug middle-end/30789] complex folding inexact

[ghazi at gcc dot gnu dot org]

------- Comment #5 from ghazi at gcc dot gnu dot org  2009-08-14 05:58 -------
Patch posted here:
http://gcc.gnu.org/ml/gcc-patches/2009-08/msg00728.html


-- 


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

[Bug middle-end/30789] complex folding inexact

[ghazi at gcc dot gnu dot org]

------- Comment #6 from ghazi at gcc dot gnu dot org  2009-08-14 16:44 -------
Subject: Bug 30789

Author: ghazi
Date: Fri Aug 14 16:44:36 2009
New Revision: 150760

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=150760
Log:
        PR middle-end/30789

        * builtins.c (do_mpc_arg2): Make extern, define for any MPC
        version.  Move declaration...
        * real.h (do_mpc_arg2): ... here.
        * fold-const.c (const_binop): Use MPC for complex MULT_EXPR
        and RDIV_EXPR.

testsuite:
        * gcc.dg/torture/builtin-math-7.c: New.


Added:
    trunk/gcc/testsuite/gcc.dg/torture/builtin-math-7.c
Modified:
    trunk/gcc/ChangeLog
    trunk/gcc/builtins.c
    trunk/gcc/fold-const.c
    trunk/gcc/real.h
    trunk/gcc/testsuite/ChangeLog


-- 


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

[Bug middle-end/30789] complex folding inexact

[ghazi at gcc dot gnu dot org]

------- Comment #7 from ghazi at gcc dot gnu dot org  2009-08-24 14:43 -------
Joseph - back in comment#2, you noted that the results of infinity and zero can
be ambiguous.  I.e. Infinity in either part of a complex number (perhaps
infinity of either sign?) is sufficient, and a pair of zeros, explicitly of
either sign is also sufficient for "zero".

This brings up the question that it's possible for GCC to produce a
compile-time result (via MPC) that is different from the runtime result (via
libgcc2) where both are C99 compliant standard conforming values.  Which one
the user receives would depend on the context (e.g. static init) or on whether
optimizations allowed GCC to fold it at compile-time.

Now this happens all the time for finite values, MPFR and MPC are more exact
that glibc's math library and often produce different results which are clearly
better.  However having GCC be more accurate than glibc (where the C library is
outside out control) is different IMHO than getting two entirely different
results from two parts of GCC, i.e. compile-time folding vs libgcc2.  E.g. (Inf
NaN) vs (0 -Inf) are both infinities per C99 Annex G.

While MPC has committed to producing C99 conforming results for the special
cases covered in Annex G, it is entirely possible they may come to different
conclusions as to what is best where the standard leaves things ambiguous.

Should I continue to pursue having GCC fold the Annex G special cases via MPC
if it leads to this kind of discrepancy?  Or should be seek to fold these
internally, essentially dulpicating the libgcc2 algorithm using real.c
functions so that both produce the same result?


-- 


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

[Bug middle-end/30789] complex folding inexact

[joseph at codesourcery dot com]

------- Comment #8 from joseph at codesourcery dot com  2009-08-24 16:59 -------
Subject: Re:  complex folding inexact

On Mon, 24 Aug 2009, ghazi at gcc dot gnu dot org wrote:

> This brings up the question that it's possible for GCC to produce a
> compile-time result (via MPC) that is different from the runtime result (via
> libgcc2) where both are C99 compliant standard conforming values.  Which one
> the user receives would depend on the context (e.g. static init) or on whether
> optimizations allowed GCC to fold it at compile-time.
> 
> Now this happens all the time for finite values, MPFR and MPC are more exact
> that glibc's math library and often produce different results which are clearly
> better.  However having GCC be more accurate than glibc (where the C library is
> outside out control) is different IMHO than getting two entirely different
> results from two parts of GCC, i.e. compile-time folding vs libgcc2.  E.g. (Inf
> NaN) vs (0 -Inf) are both infinities per C99 Annex G.

There are certainly other cases where folding may produce different 
results from those produced at runtime.  If a processor may require kernel 
support for subnormals or other special values we fold following IEEE 
rather than having options to declare the exact state of kernel support.  
We do not necessarily always produce a NaN with the same significand as 
hardware does.  Out-of-range conversions from floating-point to integer 
(unspecified value in Annex F) do not necessarily produce consistent 
results between folding and runtime (or between software and hardware 
floating point).  If we implement IBM long double constant folding (bug 
26374) I expect that will produce the nearest representable value to the 
exact result rather than always producing the same value as at runtime 
which might sometimes be further away from the exact result.

> Should I continue to pursue having GCC fold the Annex G special cases via MPC
> if it leads to this kind of discrepancy?  Or should be seek to fold these

I think folding these cases via MPC is fine.


-- 


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

[Bug middle-end/30789] complex folding inexact

[ghazi at gcc dot gnu dot org]

------- Comment #9 from ghazi at gcc dot gnu dot org  2009-08-31 03:05 -------
Patch for remaining issues posted here:
http://gcc.gnu.org/ml/gcc-patches/2009-08/msg01614.html


-- 


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

[Bug middle-end/30789] complex folding inexact

[ghazi at gcc dot gnu dot org]

------- Comment #10 from ghazi at gcc dot gnu dot org  2009-09-20 15:39 -------
Subject: Bug 30789

Author: ghazi
Date: Sun Sep 20 15:39:22 2009
New Revision: 151904

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=151904
Log:
        PR middle-end/30789
        * builtins.c (do_mpc_arg2): Accept DO_NONFINITE parameter.
        (do_mpc_ckconv): Accept FORCE_CONVERT parameter.
        (fold_builtin_2, do_mpc_arg1): Update accordingly.
        * fold-const.c (const_binop): Likewise.
        * real.h (do_mpc_arg2): Update prototype.

testsuite:
        * gcc.dg/torture/builtin-math-7.c: Update for testing Annex G
        cases in static initializers.


Modified:
    trunk/gcc/ChangeLog
    trunk/gcc/builtins.c
    trunk/gcc/fold-const.c
    trunk/gcc/real.h
    trunk/gcc/testsuite/ChangeLog
    trunk/gcc/testsuite/gcc.dg/torture/builtin-math-7.c


-- 


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

[Bug middle-end/30789] complex folding inexact

[ghazi at gcc dot gnu dot org]

------- Comment #11 from ghazi at gcc dot gnu dot org  2009-09-20 16:08 -------
Fixed, but depends on hard-requiring MPC.


-- 

ghazi at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
  BugsThisDependsOn|                            |40302
             Status|ASSIGNED                    |RESOLVED
      Known to work|                            |4.5.0
         Resolution|                            |FIXED
   Target Milestone|---                         |4.5.0


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

[Bug middle-end/30789] complex folding inexact

[ghazi at gcc dot gnu dot org]

------- Comment #12 from ghazi at gcc dot gnu dot org  2009-12-06 16:11 -------
Subject: Bug 30789

Author: ghazi
Date: Sun Dec  6 16:11:06 2009
New Revision: 155023

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=155023
Log:
        PR middle-end/30447
        PR middle-end/30789
        PR other/40302

        * configure.ac: Require MPC.
        * configure: Regenerate.
gcc:
        * doc/install.texi: Document MPC is required.


Modified:
    trunk/ChangeLog
    trunk/configure
    trunk/configure.ac
    trunk/gcc/ChangeLog
    trunk/gcc/doc/install.texi


-- 


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