Bug 24581 - Complex arithmetic on special cases is incorrect.
Summary: Complex arithmetic on special cases is incorrect.
Status: RESOLVED FIXED
Alias: None
Product: gcc
Classification: Unclassified
Component: c (show other bugs)
Version: 3.4.4
: P2 normal
Target Milestone: 4.5.0
Assignee: Not yet assigned to anyone
URL:
Keywords: wrong-code
: 43639 (view as bug list)
Depends on:
Blocks: 16989
  Show dependency treegraph
 
Reported: 2005-10-29 19:27 UTC by kargls
Modified: 2010-11-23 01:39 UTC (History)
8 users (show)

See Also:
Host:
Target:
Build:
Known to work:
Known to fail:
Last reconfirmed: 2009-01-29 21:36:23


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description kargls 2005-10-29 19:27:46 UTC
/*
 * The C99 standard intends x+I*y to be used for this, but x+I*y is
 * currently unusable because gcc introduces many overflow,
 * underflow, sign and efficiency bugs by rewriting I*y as
 * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
 * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
 * to -0.0+I*0.0.
 */
#include <complex.h>
#include <math.h>
#include <stdio.h>

int main(void) {

  double complex z;
  double x, y;

  x = 0.;
  y = 1. / x;
  x = copysign(x, -1.);

  /* z = 0 + i (-0) */
  z = I * x;
  printf("%e %e\n", creal(z), cimag(z));
 
  /* z = 0 + i Inf */
  z = I * y;
  printf("%e %e\n", creal(z), cimag(z));
 
}

kargl[223] gcc -o z a.c -lm
kargl[224] ./z
-0.000000e+00 0.000000e+00
nan inf

This bug is in 3.4.4 up to an including mainline.
Comment 1 jsm-csl@polyomino.org.uk 2005-10-29 20:09:45 UTC
Subject: Re:   New: Complex arithmetic on special cases is
 incorrect.

On Sat, 29 Oct 2005, kargl at gcc dot gnu dot org wrote:

>  * underflow, sign and efficiency bugs by rewriting I*y as
>  * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.

Note that the correct form is (0.0+I)*y, since I is (per C99+TC1+TC2) 
_Complex_I, of complex type (Annex G imaginary types conflicting with the 
normative standard, unless and until anything changes in this respect 
following DR#323).  But the usual arithmetic conversions as specified in 
the standard do not convert both operands to complex, so one can be real 
and one complex.

I suspect there are lots of presumptions in the compiler that arithmetic 
operations such as PLUS_EXPR and MULT_EXPR have both operands of the same 
type, which would need to be fixed to represent a real*complex product 
properly.

Comment 2 Andrew Pinski 2005-10-29 20:10:54 UTC
Hmm, I want to say this is a defect in the standard but what do I know.
Comment 3 Andrew Pinski 2005-10-29 20:14:41 UTC
(In reply to comment #1)
> I suspect there are lots of presumptions in the compiler that arithmetic 
> operations such as PLUS_EXPR and MULT_EXPR have both operands of the same 
> type, which would need to be fixed to represent a real*complex product 
> properly.

Or get the front-end emitting the correct code in the first place and not worry about real*complex in the middle-end which is a much better option and safer one.
Comment 4 Steve Kargl 2005-10-29 20:27:13 UTC
Subject: Re:  Complex arithmetic on special cases is incorrect.

On Sat, Oct 29, 2005 at 08:09:45PM -0000, joseph at codesourcery dot com wrote:
> 
> >  * underflow, sign and efficiency bugs by rewriting I*y as
> >  * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
> 
> Note that the correct form is (0.0+I)*y, since I is (per C99+TC1+TC2) 
> _Complex_I, of complex type (Annex G imaginary types conflicting with the 
> normative standard, unless and until anything changes in this respect 
> following DR#323).  But the usual arithmetic conversions as specified in 
> the standard do not convert both operands to complex, so one can be real 
> and one complex.
> 

If I read Annex G correctly, the z = I*inf = NaN + I inf
is going to really bad things because the NaN is going to
propagate if z is used in further computations.  Annex G
says z is an infinity.

Comment 5 Andrew Pinski 2005-10-29 20:51:09 UTC
Confirmed.
Comment 6 Joseph S. Myers 2009-05-08 10:22:29 UTC
Subject: Bug 24581

Author: jsm28
Date: Fri May  8 10:22:08 2009
New Revision: 147281

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=147281
Log:
	PR c/24581
	* c-typeck.c (build_binary_op): Handle arithmetic between one real
	and one complex operand specially.
	* tree-complex.c (some_nonzerop): Do not identify a real value as
	zero if flag_signed_zeros.

testsuite:
	* gcc.dg/torture/complex-sign.h: New header.
	* gcc.dg/torture/complex-sign-add.c,
	gcc.dg/torture/complex-sign-mixed-add.c,
	gcc.dg/torture/complex-sign-mixed-div.c,
	gcc.dg/torture/complex-sign-mixed-mul.c,
	gcc.dg/torture/complex-sign-mixed-sub.c,
	gcc.dg/torture/complex-sign-mul.c,
	gcc.dg/torture/complex-sign-sub.c: New tests.

Added:
    trunk/gcc/testsuite/gcc.dg/torture/complex-sign-add.c
    trunk/gcc/testsuite/gcc.dg/torture/complex-sign-mixed-add.c
    trunk/gcc/testsuite/gcc.dg/torture/complex-sign-mixed-div.c
    trunk/gcc/testsuite/gcc.dg/torture/complex-sign-mixed-mul.c
    trunk/gcc/testsuite/gcc.dg/torture/complex-sign-mixed-sub.c
    trunk/gcc/testsuite/gcc.dg/torture/complex-sign-mul.c
    trunk/gcc/testsuite/gcc.dg/torture/complex-sign-sub.c
    trunk/gcc/testsuite/gcc.dg/torture/complex-sign.h
Modified:
    trunk/gcc/ChangeLog
    trunk/gcc/c-typeck.c
    trunk/gcc/testsuite/ChangeLog
    trunk/gcc/tree-complex.c

Comment 7 Joseph S. Myers 2009-05-08 10:34:34 UTC
Mixed real/complex arithmetic now handles signed zeros properly and GCC
will no longer try to second-guess complex/complex arithmetic as having
one half real or imaginary just because the imaginary or real part of
that half is zero, so signed zeros should be handled correctly within
the constraints of not having imaginary types.  This may of course not
be what you want in that I is of complex type, not imaginary, but
imaginary types have ABI implications and are of very doubtful utility
apart from these corner cases.
Comment 8 H.J. Lu 2010-03-03 22:30:48 UTC
*** Bug 43251 has been marked as a duplicate of this bug. ***
Comment 9 H.J. Lu 2010-03-03 22:49:41 UTC
*** Bug 43251 has been marked as a duplicate of this bug. ***
Comment 10 Joseph S. Myers 2010-05-22 18:18:15 UTC
*** Bug 43639 has been marked as a duplicate of this bug. ***
Comment 11 marco atzeri 2010-11-21 19:48:40 UTC
the compiler produce incorrect output also when multiplying 
pure complex numbers (but not adding them). 
 
Using gcc (GCC) 4.5.1 20100924 (Red Hat 4.5.1-4) on  x86_64

The outcome of the following code is 

(inf,0)
(-nan,inf)
(inf,-nan)

instead of the expected

(inf,0)
(0,inf)
(inf,0)

---------------------------------------------------
#include <iostream>
#include <complex>
using namespace std;

int main()
{
        complex<double> z;
        complex<double> z2;
        complex<double> z3;
        double a = 0;
        double b = 1. / a;
        z = complex<double> (b,a);
        z2 = complex<double> (0,1);
        z3 = complex<double> (1,0);
        std::cout << z << '\n';
        z2 = z * z2 ;
        std::cout << z2 << '\n';
        z3 = z * z3 ;
}
Comment 12 Paolo Carlini 2010-11-21 22:15:32 UTC
Note that this specific PR is about *C* not C++. And the issue is supposed to be "RESOLVED FIXED". Thus, I would suggest first trying to reproduce the problem in C too and then either reopen this one or a C++ version (search Bugzilla first for duplicates).
Comment 13 marco atzeri 2010-11-21 23:16:45 UTC
(In reply to comment #12)
> Note that this specific PR is about *C* not C++. And the issue is supposed to
> be "RESOLVED FIXED". Thus, I would suggest first trying to reproduce the
> problem in C too and then either reopen this one or a C++ version (search
> Bugzilla first for duplicates).

Sorry Paolo,
I am a bit confused.

If the bug is "RESOLVED FIXED" why on 4.5.1 the outcome of the original
program is still

-0.000000e+00 0.000000e+00
nan inf
Comment 14 Paolo Carlini 2010-11-21 23:22:25 UTC
Yes I'm also a bit puzzled, either is just expected behavior or isn't really fixed ;) Myself I was surprised to see you just adding something to the audit trail as if it was just yet another testcase. Anyway, in the meanwhile I double checked that C does exactly the same (in the C++ front-end we have a completely similar piece of code, I'm not surprised), thus let's add in CC Joseph, and ask his opinion before re-opening.
Comment 15 jsm-csl@polyomino.org.uk 2010-11-21 23:33:48 UTC
For the original program I get

-0.000000e+00 -0.000000e+00
-nan inf

which appears correct (if one part of a complex number is an infinity, 
anything is valid for the other part and the overall value is still an 
infinity).
Comment 16 Steve Kargl 2010-11-21 23:43:10 UTC
On Sun, Nov 21, 2010 at 11:34:46PM +0000, joseph at codesourcery dot com wrote:
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=24581
> 
> --- Comment #15 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2010-11-21 23:33:48 UTC ---
> For the original program I get
> 
> -0.000000e+00 -0.000000e+00
> -nan inf
> 
> which appears correct (if one part of a complex number is an infinity, 
> anything is valid for the other part and the overall value is still an 
> infinity).
> 

The '-nan inf' is incorrect.  The correct answer is '0 inf'.
Comment 17 jsm-csl@polyomino.org.uk 2010-11-21 23:53:10 UTC
On Sun, 21 Nov 2010, sgk at troutmask dot apl.washington.edu wrote:

> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=24581
> 
> --- Comment #16 from Steve Kargl <sgk at troutmask dot apl.washington.edu> 2010-11-21 23:43:10 UTC ---
> On Sun, Nov 21, 2010 at 11:34:46PM +0000, joseph at codesourcery dot com wrote:
> > http://gcc.gnu.org/bugzilla/show_bug.cgi?id=24581
> > 
> > --- Comment #15 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2010-11-21 23:33:48 UTC ---
> > For the original program I get
> > 
> > -0.000000e+00 -0.000000e+00
> > -nan inf
> > 
> > which appears correct (if one part of a complex number is an infinity, 
> > anything is valid for the other part and the overall value is still an 
> > infinity).
> > 
> 
> The '-nan inf' is incorrect.  The correct answer is '0 inf'.

Annex G does not define the results for complex*complex multiplication to 
that level of detail, and for the complex*real multiplication we have here 
it seems entirely correct to have a NaN (sign unspecified) as the real 
part.
Comment 18 Steve Kargl 2010-11-22 00:05:02 UTC
On Sun, Nov 21, 2010 at 11:53:50PM +0000, joseph at codesourcery dot com wrote:
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=24581
> 
> Annex G does not define the results for complex*complex multiplication to 
> that level of detail, and for the complex*real multiplication we have here 
> it seems entirely correct to have a NaN (sign unspecified) as the real 
> part.
> 

We've had this discussion before.  Annex G, which I do acknowledge
as informative, states:

The * and / operators satisfy the following infinity properties for
all real, imaginary, and complex operands:319)

-- if one operand is an infinity and the other operand is a nonzero
   finite number or an infinity, then the result of the * operator
   is an infinity;

I'll note the above comes from n1124.pdf (page 468).  Perhaps,
the actual C99 standard says something else.

-nan is not an infinity.
Comment 19 jsm-csl@polyomino.org.uk 2010-11-22 00:11:39 UTC
On Mon, 22 Nov 2010, sgk at troutmask dot apl.washington.edu wrote:

> We've had this discussion before.  Annex G, which I do acknowledge
> as informative, states:
> 
> The * and / operators satisfy the following infinity properties for
> all real, imaginary, and complex operands:319)
> 
> -- if one operand is an infinity and the other operand is a nonzero
>    finite number or an infinity, then the result of the * operator
>    is an infinity;
> 
> I'll note the above comes from n1124.pdf (page 468).  Perhaps,
> the actual C99 standard says something else.
> 
> -nan is not an infinity.

That -nan is not an infinity is true but irrelevant, because "A complex or 
imaginary value with at least one infinite part is regarded as an infinity 
(even if its other part is a NaN)." (G.3), so the complex result of the 
multiplication *is* an infinity (with one part NaN and one part infinity, 
which is a valid representation of complex infinity).
Comment 20 marco atzeri 2010-11-22 14:17:24 UTC
(In reply to comment #19)
> On Mon, 22 Nov 2010, sgk at troutmask dot apl.washington.edu wrote:
> 
> 
> That -nan is not an infinity is true but irrelevant, because "A complex or 
> imaginary value with at least one infinite part is regarded as an infinity 
> (even if its other part is a NaN)." (G.3), so the complex result of the 
> multiplication *is* an infinity (with one part NaN and one part infinity, 
> which is a valid representation of complex infinity).

I guess that I was misleaded by the status FIXED.
Following your reasoning INVALID or WONTFIX are probably more accurate
STATUS as the behaviour is not a BUG but a possible implementation.

As  0 * Inf = NaN on real/double, it follows that for complex

( 0 + I ) * Inf = 0 * Inf + I * Inf = NaN + I * Inf 

however the implementation is not symmetric as
 
( 1 + I*0) * Inf = Inf + 0 * I

Of course (Inf + 0 * I) and (NaN + I * Inf) are both complex infinities,
but the lack of symmetry is inelegant ;-)

The table at C99 G.5.1-2 seems to suggest a symmetric behaviour, of course 
IMHO
Comment 21 Paolo Carlini 2010-11-22 14:24:08 UTC
If Joseph's comments are correct, and I trust him, then FIXED is the right status, because his patch actually fixed long standing serious issues vs the letter of Annex G.