(SC22WG14.1797) Complex bakeoff - example 3
Doug Gwyn (ACISD/MCSB)
uunet!ARL.MIL!gwyn
Fri Oct 13 10:19:22 PDT 1995
Out of fairness I should say that adding a "FLT_EPSILON" pass is
a kludge the need for which could be overlooked by a programmer
if he were not already aware that something special was needed to
handle the case "at" the singularity.
On the other hand, one wouldn't necessarily get the desired answer
even with signed-zero supporting complex arithmetic unless either
(a) he got lucky, or
(b) he anticipated the need to branch on either side of
the singularity and arranged the code to achieve that.
Basically, the original example was produced under condition (b).
If the loops had instead been
for ( r = -0.6; r <= 0.6; r += 0.05 )
for ( s = -1.5; s <= 1.5; s += .1 )
print_parts( g( r + I*h(s) ) );
which is a more natural way of doing business, one would not get
both plots for the singularity (even after allowing for the change
in output order). The example already anticipated that something
special was needed to deal with that situation, which is why it
broke the "natural" r-loop into two parts and carefully arranged
for the right "sign of zero" in each part; note that it did not
similarly break the s-loop into two parts, presumably because the
programmer already "knew" what to expect and thus knew that there
would be no issue of branches for the s-parameter.
The branching behavior "at" the singularity is dependent on the
definition of g, specifically the sqrt(i*w - 1) factor, which
resolves the question of which of the two logically valid branches
of the double-valued sqrt() "function" to select by examining the
infinitesimal imaginary perturbation. It happens in this example
with the CCE specification for the branch, that the plots come out
right, but there is no way I know of to assure that in other cases
the "right" choice will automatically be made. To the contrary,
even in the current example the "wrong" thing would have happened
had the programmer chosen how to "digitize" the loop parameters
differently:
for ( m = 0; m <= M; ++m ) {
for ( n = -N; n <= N; ++n )
printparts( g( m*dr + I * h(n*ds) ) );
for ( n = -N; n <= N; ++n )
printparts( g( (m-M)*dr + I * h(n*ds) ) );
}
I grant that the CCE approach *provides the tools* for a programmer
who is cognizant of these issues to tackle branching etc. more
cleanly than by using epsilon-kludges. But it is NOT automatic.
More information about the Numeric-interest
mailing list