Complex bakeoff - example 1

[Jim Thomas]

Complex bakeoff - example 1                              10/3/95      10:33 AM

Infinities, NaNs, and signed zeros are key to IEEE arithmetic's completeness
property:  operations always produce a well-defined result, which in turn
behaves predictably in subsequent operations.  For real arithmetic, this
property has proven useful for writing simpler more robust codes.

The special values naturally enter the complex domain as values for real and
imaginary parts.  Imaginary types were introduced into "Complex C Extensions"
(chapter 6 of  X3J11's Technical Report on Numerical C Extensions) as a way of
gaining consistency in the treatment of special values and thereby capturing
the completeness of IEEE real arithmetic for complex arithmetic.  No other way
has been presented.  Proponents of Cray's proposal  for complex arithmetic
without imaginary types have steadfastly maintained that special values are
not useful in the complex domain:  "NaNs and infinities give equal
information" and "+0 and -0 provide equal information".  Now that we have an
implementation of CCE we can offer some coded examples that clearly benefit
from the informational differences between infinities and NaNs, or +0 and -0,
and ask how they might be coded for Cray's proposal.  Of course it can be
done, even without any complex arithmetic support.  But does the code become
more complicated with Cray's "Complex Type Simplified" proposal, or less
efficient?  And perhaps more importantly, what is the effort to devise the
workaround?    

This example, from W. Kahan, illustrates the informational differences between
complex infinities and NaNs.  Here's a complex function implemented with a
simple continued fraction

double_complex cf(double_complex z) 
{ 
	return 4 - 3 / (z - 2 - 1 / (z - 7 + 10 / (z - 2 - 2 / (z - 3))));
}

This implementation works just fine with CCE, even at the points 1+0i, 2+0i,
3+0i, 4+0i, 5+2i, and 5-2i which cause an intermediate divide-by-zero.  CCE
succeeds here because it handles complex infinities as one would expect: 
finite / zero = infinite  and  finite/infinite = zero .  How should this
function be coded for Cray's complex proposal, where infinities are regarded
as no better than NaNs? 

Using the "equivalent" rational expression

double_complex rf(double_complex z)
{
	return (622 + z * (-751 + z * (324 + z * (-59 + z * 4)))) /
		(112 + z * (-151 + z * (72 + z * (-14 + z))));
}

might not be acceptable because its roundoff error near many points on the
real axis is about 10 to 100 times worse than the roundoff error for the
continued fraction.

Continued fractions are fairly common in mathematics (my Abramowitz and
Stegun, "Handbook of Mathematical Functions", has 16 index references for
them).  They are the most natural way to express some functions.   Although
widely avoided by programmers because of the relative slowness of division,
continued fractions offer some important advantages.  Their roundoff error is
often less than that of a corresponding rational function.  For those with
nonintegral coefficients the very determination of coefficients for a rational
function might incur critical roundoff errors. Also, continued fractions are
immune from over/underflow limitations, which become severe for higher degree
rational functions. 

-Jim