[Bug c/24581] New: Complex arithmetic on special cases is incorrect.

[Bug c/24581] New: Complex arithmetic on special cases is incorrect.

[kargl at gcc dot gnu dot org]

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


-- 
           Summary: Complex arithmetic on special cases is incorrect.
           Product: gcc
           Version: 3.4.4
            Status: UNCONFIRMED
          Severity: normal
          Priority: P2
         Component: c
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: kargl at gcc dot gnu dot org


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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[joseph at codesourcery dot com]

------- Comment #1 from joseph at codesourcery dot com  2005-10-29 20:09 -------
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.


-- 


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

[Bug middle-end/24581] Complex arithmetic on special cases is incorrect.

[pinskia at gcc dot gnu dot org]

------- Comment #2 from pinskia at gcc dot gnu dot org  2005-10-29 20:10 -------
Hmm, I want to say this is a defect in the standard but what do I know.


-- 

pinskia at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
          Component|c                           |middle-end


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

[Bug middle-end/24581] Complex arithmetic on special cases is incorrect.

[pinskia at gcc dot gnu dot org]

------- Comment #3 from pinskia at gcc dot gnu dot org  2005-10-29 20:14 -------
(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.


-- 


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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[sgk at troutmask dot apl dot washington dot edu]

------- Comment #4 from sgk at troutmask dot apl dot washington dot edu  2005-10-29 20:27 -------
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.


-- 


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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[pinskia at gcc dot gnu dot org]

------- Comment #5 from pinskia at gcc dot gnu dot org  2005-10-29 20:51 -------
Confirmed.


-- 

pinskia at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |NEW
     Ever Confirmed|0                           |1
   Last reconfirmed|0000-00-00 00:00:00         |2005-10-29 20:51:09
               date|                            |


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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[jsm28 at gcc dot gnu dot org]

------- Comment #6 from jsm28 at gcc dot gnu dot org  2009-05-08 10:22 -------
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


-- 


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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[jsm28 at gcc dot gnu dot org]

------- Comment #7 from jsm28 at gcc dot gnu dot org  2009-05-08 10:34 -------
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.


-- 

jsm28 at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
         Resolution|                            |FIXED
   Target Milestone|---                         |4.5.0


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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[hjl dot tools at gmail dot com]

------- Comment #8 from hjl dot tools at gmail dot com  2010-03-03 22:30 -------
*** Bug 43251 has been marked as a duplicate of this bug. ***


-- 

hjl dot tools at gmail dot com changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |a dot kumar at alumni dot
                   |                            |iitm dot ac dot in


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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[hjl dot tools at gmail dot com]

------- Comment #9 from hjl dot tools at gmail dot com  2010-03-03 22:49 -------
*** Bug 43251 has been marked as a duplicate of this bug. ***


-- 


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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[jsm28 at gcc dot gnu dot org]

------- Comment #10 from jsm28 at gcc dot gnu dot org  2010-05-22 18:18 -------
*** Bug 43639 has been marked as a duplicate of this bug. ***


-- 

jsm28 at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |svfuerst at gmail dot com


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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[paolo.carlini at oracle dot com]

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

--- Comment #21 from Paolo Carlini <paolo.carlini at oracle dot com> 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.

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[marco_atzeri at yahoo dot it]

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

--- Comment #20 from marco atzeri <marco_atzeri at yahoo dot it> 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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[joseph at codesourcery dot com]

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

--- Comment #19 from joseph at codesourcery dot com <joseph at codesourcery dot com> 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).

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[sgk at troutmask dot apl.washington.edu]

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

--- Comment #18 from Steve Kargl <sgk at troutmask dot apl.washington.edu> 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.

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[joseph at codesourcery dot com]

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

--- Comment #17 from joseph at codesourcery dot com <joseph at codesourcery dot com> 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.

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[sgk at troutmask dot apl.washington.edu]

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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[joseph at codesourcery dot com]

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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[paolo.carlini at oracle dot com]

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

Paolo Carlini <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jsm28 at gcc dot gnu.org

--- Comment #14 from Paolo Carlini <paolo.carlini at oracle dot com> 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.

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[marco_atzeri at yahoo dot it]

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

--- Comment #13 from marco atzeri <marco_atzeri at yahoo dot it> 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

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[paolo.carlini at oracle dot com]

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

--- Comment #12 from Paolo Carlini <paolo.carlini at oracle dot com> 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).

[Bug c/24581] Complex arithmetic on special cases is incorrect.

[marco_atzeri at yahoo dot it]

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

marco atzeri <marco_atzeri at yahoo dot it> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |marco_atzeri at yahoo dot
                   |                            |it

--- Comment #11 from marco atzeri <marco_atzeri at yahoo dot it> 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 ;
}