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
[EMAIL PROTECTED]
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.
