/****************************************************
*    Program to demonstrate the complex domain      *
*             Allroot 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 two complex roots of zē+1 = 0   *
*                                                   *
* SAMPLE RUN:                                       *
*                                                   *
* Input the initial guesses and their bounds:       *
*                                                   *
*    X0 = 4                                         *
*    Bond on X0 = 1                                 *
*    Y0 = 4                                         *
*    Bond on Y0 = 1                                 *
*                                                   *
* Convergence criterion: 1e-6                       *
* Maximum number of iterations: 30                  *
* How many roots are to be sought: 2                *
* Is the function defined in 1000 (1)               *
* or is it a series (2): 1                          *
*                                                   *
* The estimated roots are:                          *
*                                                   *
*      X = -0.000018                                *
*      Y = 1.000032                                 *
*                                                   *
*      X = 4.020035                                 *
*      Y = 4.028881                                 *
*                                                   *
* (Note: only one root is found)                    *
****************************************************/
#include <stdio.h>
#include <math.h>

#define  PI  4*atan(1)

double A[10], X[10], Y[10];
double b1,b2,e,u,v,u1,v1,u2,v2,u3,v3,x0,y00,xx,yy,x4,y4,z1,z2;
int    i,jj1,k,m,n,n2,n3;


// Rectangular to polar conversion
void RectPol(double x, double y, double *u, double *v)  {
  *u=sqrt(x*x+y*y);
  //Guard against ambiguous vector
  if (y==0.0)  y=1e-16;
  //Guard against divide by zero
  if (x==0.0)  x=1e-16;
  *v=atan(y/x);
  //Check quadrant and adjust
  if (x<0.0) *v=*v+PI;
  if (*v<0.0) *v=*v+2.0*PI;
}

//Polar to rectangular conversion
void PolRect(double u, double v, double *x, double *y)  {
  *x=u*cos(v);
  *y=u*sin(v);
}

//Polar power
void PolPower(int n, double u, double v, double *u1, double *v1) {
  int i;
  *u1=u;
  for (i=2; i<n+1; i++) *u1 *= u;
  *v1=n*v;
  *v1=*v1-(2.0*PI)*floor(*v1/2.0/PI);
}

//Rectangular complex number power
void ComplexPower(int n, double x, double y, double *x1, double *y1) {
  double u,v,u1,v1;
  //Rectangular to polar conversion
  RectPol(x,y,&u,&v);
  //Polar power
  PolPower(n,u,v,&u1,&v1);
  //Polar to rectangular conversion
  PolRect(u1,v1,x1,y1);
}

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

/**********************************************
*    Complex series evaluation subroutine     *
* ------------------------------------------- *
* The series coefficients are A(i), assumed   *
* real. The order of the polynomial is m. The *
* routine uses repeated calls to the nth      *
* power (Z^N) complex number subroutine.      *
* INPUTS:    x,y,m and A(i)                   *
* OUTPUTS:   z1(real) and z2(imaginary)       *
**********************************************/
void Series()  {
  double a1,a2,x1,y1;
  z1=A[0]; z2=0.0;
  //Store x and y
  a1=xx; a2=yy;
  for (n=1; n<m+1; n++) {
    //Recall original x and y
    xx=a1; yy=a2;
    //Call Z^N routine
    ComplexPower(n,xx,yy,&x1,&y1);
    //Form partial sum
    z1=z1+A[n]*x1; z2=z2+A[n]*y1;
  }
  //Restore x and y
  xx=a1; yy=a2;
}

//***** Supervisor subroutine *****
void Supervisor() {
  // Label: e100
  int j2,n5; double a4,u5,v5;
  n5=n; u5=u1; v5=v1;
  //Function or series ? 
  if (n3==1)  Eval(xx,yy,&u,&v);
  if (n3==2)  Series();
  if (n3==1)  goto e100;
  u=z1; v=z2;
  //Restore parameters
  n=n5; u1=u5; v1=v5;
e100: if (jj1==0) return;
  //Divide by the jj1 roots alscanfy found
  for (j2=1; j2<jj1+1; j2++) {
    u5=u;
    u=(xx-X[j2])*u+(yy-Y[j2])*v;
    v=(xx-X[j2])*v-(yy-Y[j2])*u5;
    a4=(xx-X[j2])*(xx-X[j2])+(yy-Y[j2])*(yy-Y[j2]);
    //Gard against divide by zero
    if (a4==0.0) a4=1e-7;
    v=v/a4; u=u/a4;
  } 
} //Supervisor()

/************************************************
* Mueller"s method for complex roots subroutine *
* --------------------------------------------- *
* This routine uses the parabolic fitting tech- *
* nique associated with Mueller"s method, but   *
* does it in the complex domain.                *
* --------------------------------------------- *
* INPUTS:                                       *
*   x0, y0 - The initial guess                  *
*   b1, b2 - A bound on the error in this guess *
*   e - The convergence criteria                *
*   n - The maximum number of iterations        *
* OUTPUTS:                                      *
*   x,y - The value of the complex root found   *
*   k - The number of iterations performed,     *
************************************************/
void ZMueller()  {
  //Labels: e100,e200,e300
  double aa,x1,y1,x2,y2,x3,y3;
  double a1,a2,b,c1,c2,d,d1,d2,e1,e2,xl1,xl2;
  //Start calculations
  k=1;
  x3=x0; y3=y00; x1=x3-b1; y1=y3-b2; x2=x3+b1; y2=y3+b2;
e100: d=(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);
  //Avoid divide by zero
  if (d==0.0) d=1e-7;
  xl1=(x3-x2)*(x2-x1)+(y3-y2)*(y2-y1);
  xl1=xl1/d;
  xl2=(x2-x1)*(y3-y2)+(x3-x2)*(y2-y1);
  xl2=xl2/d;
  d1=(x3-x1)*(x2-x1)+(y3-y1)*(y2-y1);
  d1=d1/d;
  d2=(x2-x1)*(y3-y1)+(x3-x1)*(y2-y1);
  d2=d2/d;
  //Get function values
  xx=x1; yy=y1; Supervisor(); u1=u; v1=v;
  xx=x2; yy=y2; Supervisor(); u2=u; v2=v;
  xx=x3; yy=y3; Supervisor(); u3=u; v3=v;
  //Calculate Mueller parameters
  e1=u1*(xl1*xl1-xl2*xl2)-2.0*v1*xl1*xl2-u2*(d1*d1-d2*d2);
  e1=e1+2.0*v2*d1*d2+u3*(xl1+d1)-v3*(xl2+d2);
  e2=2.0*xl1*xl2*u1+v1*(xl1*xl1-xl2*xl2)-2.0*d1*d2*u2-v2*(d1*d1-d2*d2);
  e2=e2+u3*(xl2+d2)+v3*(xl1+d1);
  c1=xl1*xl1*u1-xl1*xl2*v1-d1*xl1*u2+xl1*d2*v2+u3*xl1;
  c1=c1-u1*xl2*xl2-v1*xl1*xl2+u2*xl2*d2+v2*d1*xl2-v3*xl2;
  c2=xl1*xl2*u1+xl1*xl1*v1-d2*xl1*u2-xl1*d1*v2+v3*xl1;
  c2=c2+u1*xl1*xl2-v1*xl2*xl2-u2*xl2*d1+v2*d2*xl2+u3*xl2;
  b1=e1*e1-e2*e2-4.0*(u3*d1*c1-u3*d2*c2-v3*d2*c1-v3*d1*c2);
  b2=2.0*e1*e2-4.0*(u3*d2*c1+u3*d1*c2+v3*d1*c1-v3*d2*c2);
  //Guard against divide by zero
  if (b1==0.0) b1=1e-7;
  aa=atan(b2/b1); aa=aa/2.0;
  b=sqrt(sqrt(b1*b1+b2*b2));
  b1=b*cos(aa); b2=b*sin(aa);
  a1=(e1+b1)*(e1+b1)+(e2+b2)*(e2+b2);
  a2=(e1-b1)*(e1-b1)+(e2-b2)*(e2-b2);
  if (a1>a2) goto e200;
  a1=e1-b1; a2=e2-b2; goto e300;
e200: a1=e1+b1; a2=e2+b2;
e300: aa=a1*a1+a2*a2;
  xl1=a1*d1*u3-a1*d2*v3+a2*u3*d2+a2*v3*d1;
  //Guard against divide by zero
  if (aa==0.0) aa=1e-7;
  xl1=-2.0*xl1/aa;
  xl2=-d1*u3*a2+d2*v3*a2+a1*u3*d2+a1*v3*d1;
  xl2=-2.0*xl2/aa;
  //Calculate new estimate
  xx=x3+xl1*(x3-x2)-xl2*(y3-y2);
  yy=y3+xl2*(x3-x2)+xl1*(y3-y2);
  //Test for convergence
  if (fabs(xx-x0)+fabs(yy-y00)<e) return;
  //Test for number of iterations
  if (k>=n) return;
  //Continue
  k=k+1; x0=xx; y00=yy; x1=x2; y1=y2; x2=x3; y2=y3; x3=xx; y3=yy;
  goto e100;
} //Zmueller

//***** Coefficients subroutine *****
void Coefficients()  {
  m=5; A[0]=0; A[1]=24; A[2]=-50; A[3]=35; A[4]=-10; A[5]=1;
}

/***************************************************
*     General root determination subroutine        *
* ------------------------------------------------ *
* This routine attempts to calculate the several   *
* roots of a given series or function by repea-    *
* tedly using the Zmueller subroutine and removing *
* the roots already found by division.             *
* ------------------------------------------------ *
* INPUTS: X0, Y0 - The initial guess               *
*	  b1, b2 - The bounds on this guess        *
*	  e  - The convergence criteria            *
*	  n  - The maximum number of iterations    *
*              for each root                       *
*	  n2 - the number of roots to seek         *
*      	  n3 - A flag (1) for a function F(Z)      *
*                     (2) for a polynomial.        *
* ------------------------------------------------ *
* OUTPUTS: The n2 roots found X(i), Y(i)           *
* 	   The last number of iterations, k.       *
***************************************************/
void AllRoot()  {
  //Labels: e100, e200;
  k=0;
  if (n3==1) goto e100;
  if (n3!=2) return;  //error in n3
e100: jj1=0;
  //Save the initial guess
  x4=x0; y4=y00;
  //If n3=2 get the series coefficients
  if (n3==2)  Coefficients();
  //Test for completion
e200: if (jj1==n2) return;
  //Call Zmueller subroutine
  ZMueller();
  jj1++;
  X[jj1]=xx; Y[jj1]=yy; x0=x4; y00=y4;
  //Try another pass
  goto e200;
}

void main()  {
  printf("\n Write the initial guesses and their bounds:\n\n");
  printf("    X0 = "); scanf("%lf",&x0);
  printf("    Bond on X0 = "); scanf("%lf",&b1);
  printf("    Y0 = "); scanf("%lf",&y00);
  printf("    Bond on Y0 = "); scanf("%lf",&b2);
  printf("\n Convergence criterion: "); scanf("%lf",&e);
  printf(" Maximum number of iterations: "); scanf("%d",&n);
  printf(" How many roots are to be sought: "); scanf("%d",&n2);
  printf("\n Is the function defined in Eval (1)");
  printf(" or is it a series (2) "); scanf("%d",&n3);

  AllRoot();   //Call complex domain Allroot subroutine

  printf("\n\n The estimated roots are:");
  for (i=1; i<n2+1; i++) {
    printf("\n     X = %f",X[i]);
    printf("\n     Y = %f",Y[i]);
  }
  printf("\n\n The last number of iterations was: %d\n\n",k);
}

//End of file Allroot.cpp