/****************************************************
*        Program to demonstrate the complex         *
*             root counting subroutine              *
* ------------------------------------------------- *
* Reference: BASIC Scientific Subroutines, Vol. II  *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
*                                                   *
*                C++ version by J-P Moreau, Paris.  *
*                       (www.jpmoreau.fr)           *
* ------------------------------------------------- *
* Example: Find the number of complex roots of      *
*          F(z) = z^2 + 1                           *
*                                                   *
* SAMPLE RUN:                                       *
* Where is the center of the search circle (x0,y0): *
*                                                   *
*    X0 = 0                                         *
*    Y0 = 0                                         *
*                                                   *
* What is the radius of this circle: 4              *
* How many evaluation points per quadrant: 4        *
*                                                   *
* Number of complete cycles found:  2               *
* Residual:  0                                      *
*                                                   *
****************************************************/
#include <stdio.h>
#include <math.h>

#define  PI  4*atan(1)

int    i,m,nn;
int    N[80];
double a,w,x0,yy0;


//***************************************************
// Complex function(x,y) subroutine         
void Eval(double x, double y, double *u, double *v) {
  *u=x*x-y*y+1.0;
  *v=2.0*x*y;
}
//***************************************************


/************************************************
*       Complex root counting subroutine        *
* --------------------------------------------- *
* This routine calculates the number of complex *
* roots within a circle of radius w centered on *
* (x0,y0) by counting (u,v) transitions around  *
* the circumference. The input parameters are:  *
*     w      radius of the circle               *
*     x0,y0  center of the circle               *
*     m      evaluation points per quadrant     *
* The routine returns nn, the number of roots   *
* found, and a, where a<>0 indicates a failure. *
************************************************/
void RootNum()  {
  // Label e100
  int    i;
  double x,y,u,v;
  a=PI/(2.0*m);
  // Establish the N[i] array
  for (i=1; i<4*m+1; i++) {
    x=w*cos(a*(i-1))+x0;
    y=w*sin(a*(i-1))+yy0;
    Eval(x,y,&u,&v);  // Call complex function subroutine
    if (u>=0) if (v>=0) N[i]=1;
    if (u< 0) if (v>=0) N[i]=2;
    if (u< 0) if (v< 0) N[i]=3;
    if (u>=0) if (v< 0) N[i]=4;
  }
  // Count complete cycles counterclockwise
  nn=N[1]; a=0.0;
  for (i=2; i<4*m+1; i++)  {
    if (nn==N[i]) goto e100;
    if (nn!=4) if (nn==N[i]+1) a=a-1;
    if (nn==1) if (N[i]==4) a=a-1;
    if (nn==4) if (N[i]==1) a=a+1;
    if (nn+1==N[i]) a=a+1;
    nn=N[i];
e100: ; }
  // Complete circle
  if (nn!=4) if (nn==N[1]+1) a=a-1;
  if (nn==1) if (N[1]==4) a=a-1;
  if (nn==4) if (N[1]==1) a=a+1;
  if (nn+1==N[1]) a=a+1;
  a=fabs(a); nn=(int) floor(a/4);
  a=a-4.0*int(a/4);
}


void main()  {

  printf("\n Where is the center of the search circle (x0,y0): \n\n");
  printf("    X0 = "); scanf("%lf",&x0);
  printf("    Y0 = "); scanf("%lf",&yy0);
  printf("\n What is the radius of this circle: "); scanf("%lf",&w);
  printf("\n How many evaluation points per quadrant: "); 
  scanf("%d",&m);
  
  if (m>20) m=20;
  
  RootNum();   // Call RootNum subroutine

  printf("\n\n Number of complete cycles found: %d\n\n",nn);
  if (a!=0.0) printf(" Residual: %f\n",a);
  else printf(" Residual: 0\n\n");

}


// End of file Rootnum.cpp