Complex bakeoff - example 4

[Jim Thomas]

Complex bakeoff - example 4                              10/13/95     10:02 AM

In "Augmenting a Programming Language with Complex Arithmetic" (NCEG/91-039), 
W. Kahan and I presented two examples illustrating how the sign of zero 
can be used to distinguish the two sides of a branch cut to simplify the 
conceptualization and coding for problems involving branch cuts.  Bakeoff 
example 3 covered the first of these examples.  This bakeoff example 4 covers 
the second one.  In both examples, a simple, straightforward program solves
the 
problem, but the program depends on consistent treatment of the sign of zero. 

For both examples, the paper contains a plot showing the bad affects on the 
program of treating -0 like +0.  This example better suggests the potential
challenge of debugging and fixing such a program failure.

For comparison purposes, how might one code this example for a system 
which does not support the sign of zero?  Does the code become more
complicated 
or less efficient?  What is the effort to devise the workaround?  See Harry H.
Cheng's 
"Handling of Complex Numbers in the C^H Programming Language" (X3J11.1/93-020)

for a detailed analysis of this example.

Example 4

We wish to plot  z = g(w)  for
 
        g(w) =  1 + w*w + w*sqrt(w*w + 1) + log(w*w + w*sqrt(w*w + 1))
        
as  w  runs on radial straight lies from  0  in the right half-plane, 
including the imaginary axis.

This program uses the same scheme as bakeoff example 3 for generating files 
of coordinates and using MATLAB to plot.

//------------------------------------------------
//      Borda'a Mouthpiece
//      using a C++ prototype of "Complex C Extensions"
//------------------------------------------------

#include <stdio.h>
#include <complex.h>                    // Complex C Extensions
#include <fp.h>                         // for real atan 
 
// Constants:

const int       M       = 16;           // Number of radial lines = 2*M+1
const int       N       = 60;           // Number of points per radial line
const double    pi      = 4 * atan(1.0);
const double    dtheta  = pi / (2 * M); // Angle between radial lines   
const double    drho    = 1.2 / N;      // Step size    

// Functions:

double_complex g(double_complex w)
{
        return 1 + w*w + w*sqrt(w*w + 1) + log(w*w + w*sqrt(w*w + 1));
}

void print_parts(double_complex z);     // as in bakeoff example 3 

// The main program generates points  z = g(rho*exp(i*theta)) .

int main()
{
        int m, n;
        double_complex z;
        for (m = -M; m <= M; m++)
        {
                for(n = 1; n <= N; n++)
                {
                        z = g(n * drho * exp(I * m * dtheta));
                        print_parts(z);
                }
        }
        return 0;               
}