Complex bakeoff - example 4

[Jim Thomas]

Complex bakeoff - example 4.1                            10/25/95     12:37 PM

NOTE: Bakeoff example 4.1 differs from example 4 in that it maps the positive
and negative imaginary axes themselves, instead of their approximations from
n*drho*exp(I*m*dtheta) .  This more clearly plots the boundary of the flow, 
and fixes a critical deficiency in the example, pointed out by Harry Cheng.

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 how debugging 
and fixing such a program failure can be challanging.

For comparison purposes, how might one code this examples for the Cray
proposal 
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.1

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;
        
// theta = -pi/2 ...

        for(n = 1; n <= N; n++)
        {
                z = g(I * n * (-drho));
                print_parts(z);
        }
        
// -pi/2 < theta < pi/2 ...

        for (m = -M + 1; m < M; m++)
        {
                for(n = 1; n <= N; n++)
                {
                        z = g(n * drho * exp(I * m * dtheta));
                        print_parts(z);
                }
        }

// theta = pi/2 ...

        for(n = 1; n <= N; n++)
        {
                z = g(I * n * drho);
                print_parts(z);
        }

        return 0;               
}

// To print  real(z)  and  imag(z)  to two files created with unique names ...

void print_parts(double_complex z)
{
        static FILE* restream;
        static FILE* imstream;
        static int isopen = 0;
        if (!isopen) 
        {
                char refile[L_tmpnam];
                char imfile[L_tmpnam];
                tmpnam(refile);
                tmpnam(imfile);
                restream = fopen(refile, "w");
                imstream = fopen(imfile, "w");
                printf("refile = %s\n", refile);
                printf("imfile = %s\n", imfile);
                isopen = 1;
        }
        fprintf(restream, " %17.10e", real(z));
        fprintf(imstream, " %17.10e", imag(z));
}