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.
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!
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.
(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.
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.
Patch posted here: http://gcc.gnu.org/ml/gcc-patches/2009-08/msg00728.html
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
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?
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.
Patch for remaining issues posted here: http://gcc.gnu.org/ml/gcc-patches/2009-08/msg01614.html
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
Fixed, but depends on hard-requiring MPC.
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