User account creation filtered due to spam.

Bug 30789 - complex folding inexact
Summary: complex folding inexact
Status: RESOLVED FIXED
Alias: None
Product: gcc
Classification: Unclassified
Component: middle-end (show other bugs)
Version: 4.3.0
: P3 normal
Target Milestone: 4.5.0
Assignee: Kaveh Ghazi
URL:
Keywords:
Depends on: 40302
Blocks: 16989 21718
  Show dependency treegraph
 
Reported: 2007-02-13 18:31 UTC by Joseph S. Myers
Modified: 2009-09-20 16:08 UTC (History)
2 users (show)

See Also:
Host:
Target:
Build:
Known to work: 4.5.0
Known to fail:
Last reconfirmed: 2009-07-21 17:20:48


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Joseph S. Myers 2007-02-13 18:31:30 UTC
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.
Comment 1 Kaveh Ghazi 2009-07-21 17:20:48 UTC
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!

Comment 2 joseph@codesourcery.com 2009-07-26 17:32:04 UTC
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.

Comment 3 Kaveh Ghazi 2009-08-12 22:28:09 UTC
(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.
Comment 4 joseph@codesourcery.com 2009-08-13 01:25:10 UTC
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.

Comment 5 Kaveh Ghazi 2009-08-14 05:58:53 UTC
Patch posted here:
http://gcc.gnu.org/ml/gcc-patches/2009-08/msg00728.html
Comment 6 Kaveh Ghazi 2009-08-14 16:44:56 UTC
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

Comment 7 Kaveh Ghazi 2009-08-24 14:43:27 UTC
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?
Comment 8 joseph@codesourcery.com 2009-08-24 16:59:36 UTC
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.

Comment 9 Kaveh Ghazi 2009-08-31 03:05:34 UTC
Patch for remaining issues posted here:
http://gcc.gnu.org/ml/gcc-patches/2009-08/msg01614.html
Comment 10 Kaveh Ghazi 2009-09-20 15:39:35 UTC
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

Comment 11 Kaveh Ghazi 2009-09-20 16:08:03 UTC
Fixed, but depends on hard-requiring MPC.

Comment 12 Kaveh Ghazi 2009-12-06 16:11:23 UTC
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