Complex Numbers

Complex Numbers

[Andreas Klein]

Hello

I have looked at the implementation of complex arithmetic in gcc.

I see tree problems with naive formulas for multiplication and
division that are currently in use.

1. Problems with special values. For example (infinity+ i*NaN)^2
   should be (infinity + i*0) according to IEC 60559 and not
   (NaN + i*NaN). Please see the example code in Annex G.

2. If a is large but representable with a^2=infinity. We get
   (a+ i*a)^2=(NaN+i*infinity) with the naive formula but we want
   (0+i*infinity)

3. The are cancellation problems.
   For example if we use 3 digit decimal floating point arithmetic with
   a=4.02, b=4.00, c=4.02 and d=4.00 the exact result for the
   real part would be $ac-bd = 0.1604$. But with the naive formula
   we get a * c = 16.2, b * c = 16.0 and thus a*c - b* d = 0.2.
   This make a 25% error.

The IEC 60559 standard deal only with problem 1 but not with the
problems 2 and 3. According to Annex G the "specifications (of IEC
60559) are not normative" and therefore I think we should do better
and use code that avoid all problems mentioned above. Perhaps we should
provide the IEC 60599 code through the __STDC_IEC_559_COMPLEX__ pragma
as suggested in Annex G.

I have a solution to the problems describe above and I am willing (and
hopefully able) to write a fix. How ever the resulting code would be
about 2-3 times slower than the naive code. For example on a computer
that has an fma machine instruction

p=-b*d;
r=fma(b,d,p);
x=fma(a,c,p)-r;

computes the real part avoiding cancellation, but this needs 4
operations where the naive formula needs only 2.

p=-b*d;
x=fma(a,c,p)

So please tell me do you think we should switch to more accurate
algorithms or would you prefer the old implementation?

With best regards
     Andreas Klein

    Institute for Mathematics and Computer Science
    klein@mathematik.uni-kassel.de
    University of Kassel
    Germany

PS: I am better in math than in email. So please be patient if have
    chosen the wrong format or the wrong group.

Re: Complex Numbers

[Paolo Carlini]

Paolo Carlini wrote:

> We are already aware of this issue, since you have already reported
> it ;) The relevant PR is middle-end/18902.

Forgot to add: for other issues, related in particular to multiplication,
not only division, please file appropriate Bugzilla PRs.

Thanks!
Paolo.

Re: Complex Numbers

[Paolo Carlini]

> I have looked at the implementation of complex arithmetic in gcc.

We are already aware of this issue, since you have already reported
it ;) The relevant PR is middle-end/18902.

Indeed, our plan involves enabling the (*already available*) algorithm
due to Smith. There are still some open issues, however (see "Depends on"
and "blocks" fields, for further details).

As usual, patches are *welcome*, if prepared according to our (usual)
rules and the necessary copyright assignment is in place.

Paolo.

Re: complex numbers

[Andreas Klein]

Hello
> >
> Ok, thanks. The important section of the C99 standard is Annex G (IEC
> 60559-compatible complex arithmetic): it even provides a reference
> implementation of the division in Example2. Perhaps, you could have a
> look to a public draft of the final standard, just Google a bit... ;)
> But really you should get hold of the real document, is now published by
> Wiley and most scientific libraries should have a copy of it.

All I could find is the the CX_LIMITED_RANGE pragrama can be used to allow
the compiler a bad but fast complex implementation. But I am not able to
find a copy of the C99 standard or the IEC 60559 standard. Can you send me
the Example2? Or does any one know a german (perfeable Hessen) university
where I can get these standards?

Yours
	Andreas Klein

Re: complex numbers

[Paolo Carlini]

Andreas Klein wrote:

>Unfortunally I have own no copy of the C99 standard. So I would be glade
>if you could give me an internet ressource which disscuss the C99 division
>algorithm or something like that. Then I will try to check what we can do.
>  
>
Ok, thanks. The important section of the C99 standard is Annex G (IEC 
60559-compatible complex arithmetic): it even provides a reference 
implementation of the division in Example2. Perhaps, you could have a 
look to a public draft of the final standard, just Google a bit... ;) 
But really you should get hold of the real document, is now published by 
Wiley and most scientific libraries should have a copy of it.

Thanks,
Paolo.

Re: complex numbers

[Andreas Klein]

Hello
>
> | What was the critice you mentioned above? I can not imagine a sitation in
> | which I would need the naive implementation.
>
> Oh, I got repeated complaints from users that the correct method of
> computation was slow -- look at the bugzilla archive.  I believe there
> might alos be discussions on GCC main mailing list.

It seam that i have trouble to understand how bugzilla works. It would
help me if I know the exact bug number.

> The grammar school algorithms is needed for the primary template
> because it is the one that gets instantiated, for say, complex<int>.
> Notice that such datatypes are required by LIA-3 (in the last phase of
> development).

Yes that right. I was speaking of complex<floting-point> types. For
complex<integer/rational> the grammar school algorithmen is the right one.

Yours
	Andreas Klein

Re: complex numbers

[Andreas Klein]

Hello
> >
> Have a look to the implementation: it looks like that even if we switch
> to the
> better algorithm, still we don't get fully right C99. Of course this
> last point
> must be better investigate (I'm not a floating point expert) but I
> expect someone
> replying: "let's implement C99 division right instead".
>
Unfortunally I have own no copy of the C99 standard. So I would be glade
if you could give me an internet ressource which disscuss the C99 division
algorithm or something like that. Then I will try to check what we can do.
Meanwhile I think to switch the compiler flag to get the improved
algorithmen is a step in the right direction.

Yours
	Andreas Klein

Re: complex numbers

[Gabriel Dos Reis]

Andreas Klein <klein@mathematik.uni-kassel.de> writes:

| This look like a good deal. However for floting point computations I
| prevere good results over fast results.

You're in the minority (including me :-)).

-- Gaby

Re: complex numbers

[Gabriel Dos Reis]

Andreas Klein <klein@mathematik.uni-kassel.de> writes:

| On Mon, 6 Dec 2004, Gabriel Dos Reis wrote:
| 
| >
| > As a matter of fact, the implementation of <complex> is criticized,
| > once in a while, because it does NOT use the grammar school rule you
| > present above.  However, for float, double, long double it specializes
| > to __complex__ T which is  what the compiler uses to implement complex
| > numbers in C99 and Fortran.  So, the problem is a compiler problem not
| > libstdc++ problem.
| 
| I have testet my program with a 2.95 and several 3.x versionsions of gcc.
| The lastet version I have is 3.3.3.
| All versions gave the same result.

Yup.  It is inside the compiler

| What was the critice you mentioned above? I can not imagine a sitation in
| which I would need the naive implementation.

Oh, I got repeated complaints from users that the correct method of
computation was slow -- look at the bugzilla archive.  I believe there
might alos be discussions on GCC main mailing list.

The grammar school algorithms is needed for the primary template
because it is the one that gets instantiated, for say, complex<int>.
Notice that such datatypes are required by LIA-3 (in the last phase of
development). 

We don't have overload on "concepts" so, for the moment that is what
is used for anything other thand float, double, and long double.

| But I think it can not be
| good that if I get different result for complex and real arithmetic. If
| the number are the same real numbers.
| 
| >
| > Did you look at the actual implementation?
| >
| 
| As you mentinon it if have missed the specilization at the end of
| std_complex.h. Sorry. I still think that we should have and other
| implementation for complex<floating_point>, but I cannot change the code
| of __complex__ T in the complier.

I don't think we need another implementation of
complex<floating_point>.  What we need is to improve the
compiler. Convince people that we have correctness, then speed.

-- Gaby

Re: complex numbers

[Paolo Carlini]

Andreas Klein wrote:

>>... but notice that this issue is tricky: there are computational issues
>>(we are adding
>>at least a branch for each division) and correctness issues (what about
>>C99?)
>>    
>>
>As I see it the naive formula needs
>
>6 multipications, 2 divisions and 3 additions/subtractions
>
>and the improved formua needs
>
>1 comparision, 3 multipactions, 3 divisions and 3 additions/subtractions
>
>This look like a good deal. However for floting point computations I
>prevere good results over fast results.
>  
>
Thanks for the summary: I just sent a message to the gcc list, let's see 
what comes
out... Maybe you may want to reply something there.

>What do you meen with "correctness issues"?
>  
>
Have a look to the implementation: it looks like that even if we switch 
to the
better algorithm, still we don't get fully right C99. Of course this 
last point
must be better investigate (I'm not a floating point expert) but I 
expect someone
replying: "let's implement C99 division right instead".

Paolo.

Re: complex numbers

[Andreas Klein]

Hello

> >However I think if flag_complex_divide_method = 1 fix the problem it would
> >be a good idea to set it by default.
> >
> >
> ... but notice that this issue is tricky: there are computational issues
> (we are adding
> at least a branch for each division) and correctness issues (what about
> C99?)

As I see it the naive formula needs

6 multipications, 2 divisions and 3 additions/subtractions

and the improved formua needs

1 comparision, 3 multipactions, 3 divisions and 3 additions/subtractions

This look like a good deal. However for floting point computations I
prevere good results over fast results.


What do you meen with "correctness issues"?

Andreas Klein

Re: complex numbers

[Paolo Carlini]

Paolo Carlini wrote:

> I will try to do the same as soon as possible...

I can confirm that setting flag_complex_divide_method = 1 leads to (0, 0).

Paolo.

Re: complex numbers

[Paolo Carlini]

Andreas Klein wrote:

>>Have a look to expand_complex_division in gcc/tree-complex.c, then
>>gcc/toplev.c for flag_complex_divide_method.
>>
>>Andreas, just for curiosity, are you willing to rebuild your gcc
>>with flag_complex_divide_method = 1 and report???
>>    
>>
>Willing is not the problem. But I have only limitet privileges on this
>machine and I am not sure if I can find an other machine with enough
>resources.
>  
>
I will try to do the same as soon as possible...

>However I think if flag_complex_divide_method = 1 fix the problem it would
>be a good idea to set it by default.
>  
>
... but notice that this issue is tricky: there are computational issues 
(we are adding
at least a branch for each division) and correctness issues (what about 
C99?)

Exerimenting a bit with those options and functions is only a *very* 
preliminary,
explorative step!

Paolo.

P.S. In gcc3.3.x the relevant functions are in gcc/optabs.c

Re: complex numbers

[Andreas Klein]

Hello

> Interestingly, it looks like the discussed improved algorithm is
> *already* implemented, just not used!

Curious

> Have a look to expand_complex_division in gcc/tree-complex.c, then
> gcc/toplev.c for flag_complex_divide_method.
>
> Andreas, just for curiosity, are you willing to rebuild your gcc
> with flag_complex_divide_method = 1 and report???

Willing is not the problem. But I have only limitet privileges on this
machine and I am not sure if I can find an other machine with enough
resources.

However I think if flag_complex_divide_method = 1 fix the problem it would
be a good idea to set it by default.

Andreas Klein

Re: complex numbers

[Paolo Carlini]

> As you mentinon it if have missed the specilization at the end of
> std_complex.h. Sorry. I still think that we should have and other
> implementation for complex<floating_point>, but I cannot change the code
> of __complex__ T in the complier.

Interestingly, it looks like the discussed improved algorithm is
*already* implemented, just not used!

Have a look to expand_complex_division in gcc/tree-complex.c, then
gcc/toplev.c for flag_complex_divide_method.

Andreas, just for curiosity, are you willing to rebuild your gcc
with flag_complex_divide_method = 1 and report???

Thanks,
Paolo.

Re: complex numbers

[Andreas Klein]

On Mon, 6 Dec 2004, Gabriel Dos Reis wrote:

>
> As a matter of fact, the implementation of <complex> is criticized,
> once in a while, because it does NOT use the grammar school rule you
> present above.  However, for float, double, long double it specializes
> to __complex__ T which is  what the compiler uses to implement complex
> numbers in C99 and Fortran.  So, the problem is a compiler problem not
> libstdc++ problem.

I have testet my program with a 2.95 and several 3.x versionsions of gcc.
The lastet version I have is 3.3.3.
All versions gave the same result.

What was the critice you mentioned above? I can not imagine a sitation in
which I would need the naive implementation. But I think it can not be
good that if I get different result for complex and real arithmetic. If
the number are the same real numbers.

>
> Did you look at the actual implementation?
>

As you mentinon it if have missed the specilization at the end of
std_complex.h. Sorry. I still think that we should have and other
implementation for complex<floating_point>, but I cannot change the code
of __complex__ T in the complier.

Yours
	Andreas Klein

Re: complex numbers

[Gabriel Dos Reis]

Andreas Klein <klein@mathematik.uni-kassel.de> writes:

| Hello
| 
| I have found a bug in the implementation of the complex library of
| g++ and the complex.h library of the gcc (c99 support).
| 
| The simplest program that demonstrates the bug is:
| 
| --------------------------------
| 
| #include<iostream>
| #include<complex>
| 
| using namespace std;
| 
| int main() {
|   complex<double> a = 1.0/0.0;  // + infinity
|   complex<double> b = 1.0;
|   complex<double> c = b/a;      // should be 0 but is (NaN,NaN)
| 
|   cout << a << endl;
|   cout << b << endl;
|   cout << c << endl;
| }
| 
| --------------------------------
| 
| Since all values are real one should expect the result 0
| (the IEEE floating-point value for 1/infinity). More serious is that
| if a contains a large but representable value for which |a|^2 is not
| representable the implementation will yields bad results. For example
| let a be the largest representable number and b be a/2 then b/a should
| be 0.5 and not NaN.
| 
| The bug comes from the naive implementation
| 
| a + b I    ac + bd   bc - ad
| -------  = ------- + ------- I
| c + d I    c^2+d^2   c^2+d^2
| 
| We get unnecessary overflow errors if c^2+d^2 is too large.
| 

As a matter of fact, the implementation of <complex> is criticized,
once in a while, because it does NOT use the grammar school rule you
present above.  However, for float, double, long double it specializes
to __complex__ T which is  what the compiler uses to implement complex
numbers in C99 and Fortran.  So, the problem is a compiler problem not
libstdc++ problem.

[...]

| So my suggestion is to add template specialization for the
| operator/= for the types complex<float> complex<double> and
| complex<long double> that use Smith formula. (For integer types I
| think we should stay with the old implementation.)

Did you look at the actual implementation?

-- Gaby

complex numbers

[Andreas Klein]

Hello

I have found a bug in the implementation of the complex library of
g++ and the complex.h library of the gcc (c99 support).

The simplest program that demonstrates the bug is:

--------------------------------

#include<iostream>
#include<complex>

using namespace std;

int main() {
  complex<double> a = 1.0/0.0;  // + infinity
  complex<double> b = 1.0;
  complex<double> c = b/a;      // should be 0 but is (NaN,NaN)

  cout << a << endl;
  cout << b << endl;
  cout << c << endl;
}

--------------------------------

Since all values are real one should expect the result 0
(the IEEE floating-point value for 1/infinity). More serious is that
if a contains a large but representable value for which |a|^2 is not
representable the implementation will yields bad results. For example
let a be the largest representable number and b be a/2 then b/a should
be 0.5 and not NaN.

The bug comes from the naive implementation

a + b I    ac + bd   bc - ad
-------  = ------- + ------- I
c + d I    c^2+d^2   c^2+d^2

We get unnecessary overflow errors if c^2+d^2 is too large.

For a better result one should use the following formula due to Smith
(see "What Every Computer Scientist Should Know About Floating-Point
Arithmetic"  http://cch.loria.fr/documentation/IEEE754/)

           |  a + b(d/c)   c - a(d/c)
           |  ---------- + ---------- I     if |d| < |c|
a + b I    |  c + d(d/c)   a + d(d/c)
-------  = |
c + d I    |  b + a(c/d)   -a+ b(c/d)
           |  ---------- + ---------- I     if |d| >= |c|
           |  d + c(c/d)   d + c(c/d)

Even better we should test for |d| = |c| and definie

a + b I       a + b s   c - a s
-------  =    ------- + -------   if |d| = |c|
c + d I       c + d s   d + c s

where s = 1 if c and d have the same sign bit and s = -1
otherwise. (Especially if c = +0.0 and d = -0.0 we must set s = -1.)
This improves the results if |c|=|d|=0 or |c|=|d|=infinity.

So my suggestion is to add template specialization for the
operator/= for the types complex<float> complex<double> and
complex<long double> that use Smith formula. (For integer types I
think we should stay with the old implementation.)

I would volunteer to fix the c++ library.
But I have less experience in C99 so it would be better if another can
fix the c library. I guess that the same bug will be present in the
libraries of the other languages (Fortan, etc.), but I have not tested it.

I am sorry but I have difficulties with the bug database, so I send
you the mail directly.

With best regards
     Andreas Klein

    Institute for Mathematics and Computer Science
    klein@mathematik.uni-kassel.de
    University of Kassel
    Germany