/**************************************************** * 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 #include #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