Complex bakeoff - example 3

[Harry H. Cheng]

> Complex bakeoff - example 3                              10/12/95     3:27 PM
> 
> Bakeoff example 1 used complex infinity.  Example 2 used multiple 
> complex infinities.  This example uses the sign of zero.
> 
> 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.  Below
> is a simple, straightforward C program for the first of the two examples, 
> but the program depends on consistent treatment of the sign of zero.  The
> paper includes plots showing adverse affects of treating -0 like +0.  
> 
> This example is intended to illustrate a basic use of the sign of zero.  More 
> complicated examples better dramatize the relative difficulty in developing 
> a solution that does not utilize the sign of zero.  Still, 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?

I have tested your program with my implementation with sign
of zeros being ignored in complex number.
But, when I test the output, I use my plot routines.

Besides the change from double_complex to double complex,
The workaround I made to your program is noted in your
following original program.
I do not think the code is much more complicated.

As pointed out in my paper,
it appears that the singularity of the slit along the imaginary axis between
1<z^2<0 is the source of the problem, the images for the segments 
|y|>1 along the imaginary axis in the z-plane are well behaved.

 
> 
> In "Handling of Complex Numbers in the C^H Programming Language" 
> (X3J11.1/93-020), section 7.7, Harry H. Cheng addresses these very examples, 
> with reference to his experimental extended C implementation.  He rightly 
> concludes that distinguishing the signs of zero is not necessary to program 
> the two examples.  However, his schemes are more complicated, and obviously 
> required time to devise and prove.
> 
> Example 3
> 
> As a conformal map,  w = f(z) = (z - 1/z) / 2  maps the complex z-plane 
> twice, once for  |z| >= 1  and once for   |z| <= 1 , onto the complex w-plane,
> 
> mapping the unit circle  |z| = 1  to a slit along the imaginary axis from  
> w = -i  to w = +i.  The inverse map
> 
>         g(w) = w - i*sqrt(i*w - 1)*sqrt(i*w + 1)
> 
> maps the whole w-plane slitted along  w = -i  to w = +i  onto the outside of
> the 
> z-plane's unit circle  |z| >= 1 .  Vertical lines in the w-plane map to the 
> streamlines of a vertical "eluding" flow around the unit circle on the
> z-plane.
> We wish to calculate points to exhibit these streamlines.
> 
> To smooth the plot near stagnation points  g(+/-i)  where  g'(+/-i)  is
> infinite
> we use
> 
>         h(s) = s*(s*s*(3*s*s - 10) + 15) / 8
> 
> which satisfies  h(-s) = -h(s) ,  h'(s) >= 0 ,  h(1) = 1 ,  and h'(1) = 0 .
> For any fixed real  r ,  w = r + i*h(s)  traces out a vertical straight line 
> segment as we let  s  run through  -1.5 <= s <= 1.5  , in small steps
> ds = 3/32.  z = g(w)  traces out a streamline past the circle.
>         
> To plot several streamlines, we run  r  through  0 <= r <= 0.6 , in steps of 
> dr = 0.05 .  For each such  r  we plot two streamlines,  z = g(r + i*h(s))
> to the right of the circle and  z = g(-r + i*h(s))  to the left.
> 
> The main program calculates the points  z  on the streamlines and prints 
> the real and imaginary parts into two files.  MATLAB can be used to load 
> the file of imaginary parts and the file of real parts into array variables, 
> reshape them into  (2*N+1)x(2*M+2)  arrays  X  and  Y , respectively, and 
> then plot column of  Y   against column of  X .
> 
> //-------------------------------------------------------
> //      Eluding flow past a disk
> //      using a C++ prototype of "Complex C Extensions"
> //-------------------------------------------------------
> 

  #include <float.h>

> #include <stdio.h>
> #include <complex.h>                    // Complex C Extensions
> 
> // Constants:
> 
> const int       M       = 12;           // Number of streamlines = 2*M+1 (one
> splits) 
> const int       N       = 16;           // Number of points per steamline =
> 2*N+1
> const double    dr      = 0.05;         // Distance between streamlines
> const double    ds      = 3.0 / 32.0;   // Step-size

  double complex I = complex(0,1);

> 
> // Functions:
> 
> double_complex g(double_complex w) { return w - I * sqrt(I * w - 1) * sqrt(I *
> w + 1); }
> double h(double s) { return s * (s * s * (3 * s * s - 10) + 15) / 8 ; }
> void print_parts(double_complex z);
> 
> // To generate points along streamlines ...
> 
> void main()
> {
>         double_complex z;
>         int m,n;
>         for (m = 0; m <= M; m++) 
>         {
>                 // right-side streamline 
>                 for (n = -N; n <= N; n++) 
>                 {
>                         z = g(m * dr + I * h(n * ds));
>                         print_parts(z);
>                 }               
>                 // left-side streamline 
>                 for (n = -N; n <= N; n++) 
>                 {
>                         z = g(-(m * dr) + I * h(n * ds));
>                         print_parts(z);
>                 }    
           
                  // get the missing left slit of singularity points
                  // which appear on right-side
                  for (n = -N; n <= N; n++) 
                  {
                         z = g(-FLT_EPSILON + I * h(n * ds));
                         print_parts(z);
                  } 
              
>         }
> }
> 
> // 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));
> }